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COMPUTER PROGRAM FOR DESIGN OF TWO-DIMENSIONAL SHARP-EDGED -THROAT 

SUPERSONIC NOZZLE WITH BOUNDARY -LAYER CORRECTION 
by Louis J. Goldman and Michael R. Vanco 
Lewis Research Center 

SUMMARY 

A FORTRAN IV computer program is presented for the design of sharp-edged- 
throat supersonic nozzles where losses are accounted for by correcting the ideal nozzle 
geometry for boimdary- layer displacement thickness. The ideal nozzle is designed by 
the method of characteristics. Boundary- layer parameters are calculated by Cohen and 
Reshotko’s method for laminar flow and by Sasman and Cresci’s method for turbulent 
flow. The final nozzle geometry is then obtained by adding the displacement thicknesses 
to the ideal nozzle coordinates. The program is essentially a modified combination of 
two programs presented previously. 

The computer program input consists primarily of the exit Mach number, specific 
heat ratio, and total flow conditions. The program output gives the corrected nozzle 
geometry. An example of the output is included to indicate the use of the program. 


INTRODUCTION 

Methods for the design of supersonic turbines for possible use in turbopump and 
open-cycle auxiliary power systems where high-energy fluids are used and high pres- 
sure ratios are available have recently become of interest. The method of character- 
istics as applied to the two-dimensional isentropic flow of a perfect gas can be used for 
the design of the supersonic blading. Computer programs for the isentropic design of 
two-dimensional supersonic nozzles and rotor blade sections have been described in 
reference 1 and 2, respectively. The design of blading by these procedures must then 
be supplemented by a knowledge of the losses that occur in the nozzle and rotor. 

The purpose of this report is to present a computer program for the design of 
sharp-edged-throat supersonic nozzles where losses are accounted for by correcting 
the ideal nozzle geometry for boundary- layer displacement thickness. The ideal nozzle 



geometry is obtained using the computer program described in reference 1. Boundary- 
layer parameters are calculated using the computer program described in reference 3. 
The final nozzle geometry is then obtained by adding the displacement thicknesses to the 
ideal nozzle coordinates. The program described herein is essentially a modified com- 
bination of the two programs presented in references 1 and 3 . 

The boxmdary- layer parameters are also used to calculate the conditions downstream 
of the nozzle after the flow has mixed to a uniform state. The procedure described in 
reference 4 is used for this purpose. 

In this report, a description of the input and output and a complete FORTRAN IV 
listing of the program are presented. A brief description of the computer program and 
method of design are also given. An example of the output is included to indicate the 
use of the program. 


METHOD OF ANALYSIS 

The design of sharp-edged-throat supersonic nozzles where losses are accounted 
for by correcting the ideal nozzle geometry for boundary- layer displacement thickness 
is described herein. This type of ideal nozzle is desirable since it produces uniform 
parallel flow at the exit in the smallest possible distance. The ideal symmetric nozzle 
is designed by the method of characteristics as applied to the isentropic flow of a per- 
fect gas. The boundary -layer parameters (displacement and momentum thicknesses) 
are then calculated for the ideal nozzles. The final nozzle geometry is obtained by 
adding the displacement thicknesses to the ideal nozzle coordinates. In addition, the 
boundary -layer parameters are used to calculate the aftermixing conditions downstream 
of the nozzle assuming the flow mixes to a uniform state. 


Nozzle Description 

A typical supersonic turbine nozzle, as seen in figure 1, consists of three sections; 
(1) a converging (subsonic) section, (2) a diverging (supersonic) section, and (3) a 
straight section on the suction surface. The converging section produces the flow turn- 
ing with little losses and is not designed by the program. The diverging section accel- 
erates the flow to the desired free-stream Mach number at the exit. This section is 
designed by the method of characteristics using the program described in reference 1. 
The straight section on the suction surface completes the nozzle profile, and its length 
is determined by the required nozzle angle (see fig. 1). 
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Ideal nozzle 



3 



Boundary “Layer Calculations 


The boundary- layer parameters (displacement and momentum thicknesses) are 
calculated for the ideal nozzle using the computer program described in reference 3. 

The program uses Cohen and Reshotko's method (ref. 5) for the calculation of laminar 
boundary layers and Sasman and Cresci's method (ref. 6) for turbulent boundary layers. 
Curvature effects are not considered in these calculations. 

In the laminar regime, a single ordinary differential equation (the momentum inte- 
gral equation) is solved numerically. As explained in reference 5, the results have to 
be extended for flows in highly favorable pressure gradients as would occur in a nozzle. 
For turbulent flow, coupled first-order ordinary differential equations (the momentum 
and moment-of-momentum integral equations) are solved using Runge-Kutta techniques. 
Transition from laminar to turbulent flow, if it occurs, is predicted by the program. 

The displacement thicknesses obtained from the program are then added to the ideal 
nozzle coordinates to obtain the final nozzle profile. Figure 1 shows a nozzle designed 
in this manner. The dashed line represents the ideal nozzle profile. 


Aftermixing Conditions 

The boundary- layer parameters (displacement and momentum thickness) at the noz- 
zle exit are also used to calculate the aftermixing conditions by the method given by 
Stewart in reference 4. In this loss model, the flow sufficiently downstream of the blade 
row is assumed to be mixed to a uniform condition. Application of the continuity, 
momentum, and energy equations between stations 1 and 2 (fig. 1) results in the deter- 
mination of the aftermixing conditions. 

Subsonic and supersonic aftermixing axial Mach number solutions are possible for 
this loss model when the free-stream axial Mach number at the nozzle exit (before mix- 
ing) is supersonic. The subsonic solution corresponds to mixing plus oblique shock 
lossep, whereas the supersonic solution corresponds to shockless mixing. A more de- 
tailed discussion of the different solutions can be found in reference 7. 

DESCRIPTION OF INPUT 

A description of the input for the FORTRAN IV computer program is given in this 
section. The input quantities consist essentially of the nozzle-exit Mach number, 
specific-heat ratio, nozzle angle, throat half-height, nozzle subsonic section coordin- 
ates and corresponding pressure ratios, total temperatime, total pressure, gas constant, 
and initial momentum or displacement thickness. The program gas properties are set 
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up for air. For gases other than air the changes required to the program are as de- 
scribed in appendix A. 

The input format is shown in table I. The input variables are the following: 


ME 

DV 

GAM 

NP 


NS 

ALPl 

YTH 

XCP 

YCP 

PSPT 

R 

PTZ 

TTZ 

UPMACH 

NVP 

NTURB 


nozzle-exit Mach number 

increment in flow turning (recommended value, 0. 1), deg 
specific-heat ratio 

nozzle coordinate printout indicator for supersonic section (e.g. , if NP = 
10, every tenth coordinate is printed out; it should also be noted that for 
any value of NP, the last coordinate is always printed out, and NP must 
be an integer and right adjusted) 

number of input points for nozzle inlet including the throat (must have at 
least one point, the throat point) 

nozzle angle measured from tangential direction, deg, see fig. 1 
throat half-height, m (ft) 

array of X-coordinates of nozzle subsonic section (based on throat half- 
height of unity) 

array of Y-coordinates of nozzle subsonic section (based on throat half- 
height of unity) 

array of static to total pressure ratios corresponding to the nozzle subsonic 
section coordinates XCP and YCP 

gas constant, J/(kg)(K) ((ft)(lbf)/(slug)(°R)) 

p p 

inlet or upstream total pressure, N/m (Ibf/ft ) 
inlet or upstream total temperatxire, K (°R) 
inlet or upstream Mach number 

integer number of points desired in velocity profile of boundary layer at 
each station (must have at least 1) 

integer number of station, if any, at which user wishes turbulent boundary 
layer to begin (NTURB is usually zero, allowing the program to calculate 
position of transition to turbulent boundary layer; if NTURB = 1, initial 
values must be given for DTURB and TTURB) 
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KEM 

DLAM 

TLAM 

DTURB ‘ 

TTURB 

KPRE 

KGRAD 

KSDE 

KLAM 

KM AIN 

KPROF 

TE 


integer (0 or 1) indicating which of two allowable sets of units are used in 


input; 

U.S. Customary (Ibf, slugs, ft, sec, °R, and ft- lb) 0 

International System (N, kg, m, sec, K, and J) 1 


initial displacement thickness, if any, of laminar boundary layer at first 
point, m (ft) (DLAM may be zero, the stagnation point, or have some 
finite value) 

initial momentum thickness, if any, of laminar boundary layer at first point, 
m (ft) (TLAM may be zero, or have some finite value) 

initial displacement thickness, if any, of turbulent boundary layer, m (ft) 

initial momentum thickness, if any, of turbulent boundary layer, m (ft) 

integer (0 or 1) indicating whether printing of output from PRECAL is 


desired: 

Output suppressed 0 

Output printed 1 


integer (0 or 1, see KPRE) indicating whether printing of surface gradients 
of velocity and Mach number is desired 

integer (0 or 1, see KPRE) indicating whether printing of solutions of lami- 
nar and turbulent differential equations is desired 

integer (0 or 1, see KPRE) indicating whether printing of laminar calcula- 
tions for location of instability and transition is desired 

integer (0 or 1, see KPRE) indicating whether printing of principal calcula- 
ted boundary- layer parameters is desired 

integer (0 or 1, see KPRE) indicating whether printing of velocity profiles 
is desired 

thickness of trailing edge, m (ft) 
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TABLE I. - INPUT FORM 


[Numbers in corners are card column numbers. ] 


1 6 
ME 

7 12 

DV 

13 18 

GAM 

20 

NP 

22 

NS 

23 28 

ALPl 

29 38 

YTH 

1 




1 

1 



1 

10 

XCP 

20 

YCP 

30 

PSPT 

40 





10 

R 

20 

PTZ 

30 

TTZ 

40 

UP MACH 





5 

NVP 

10 

NTURB 

15 

KEM 





10 

DLAM 

20 

TLAM 

30 

DTURB 

40 

TTURB 





5 

KPRE 

10 

KGRAD 

15 

KSDE 

20 

KLAM 

25 

KMAIN 

30 

KPROF 









6 

TE 





SPECIAL INSTRUCTIONS FOR PREPARING INPUT 
Subsonic Section 

For this program, the subsonic section of the nozzle is used to determine the 
boundary- layer parameters at the throat. If the displacement thickness or momentiun 
thickness at the throat is known, it is put in as input, DLAM or TLAM. Then, only one 
point is required for the subsonic section, the throat point. At this point, X = 0, Y = 1.0, 
and PSPT equals the static to total pressure ratio at Mach 1. 

If the boundary- layer parameters are not known, a subsonic section is put in. Fig- 
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Figure 2. - Coordinate system for geometry input. Ali coordinates referenced to throat 
haif-height of unity. 


ure 2 shows an example subsonic section with the required coordinate system. The sub- 
sonic section must be symmetrical because only one surface is used, and the X-Y 
coordinates must be referenced to a throat half-height of unity. Also, along with coor- 
dinates, the static to total pressure ratio (PSPT) at each point is needed. 


Output 

Usually KPRE, KGRAD, KSDE, KLAM, and KPROF are set equal to zero. They 
do not give the main output of the boundary- layer section. However, if this additional 
output is desired, these quantities are set equal to 1. A description of this output is 
given in appendix B. 


DESCRIPTION OF AAAIN OUTPUT 

An example of the output from the program is shown in table II . This output is in 
U.S. Customary units. Each section of the output has been numbered to correspond to 
the following description: 

(1) The first output of the program is the input values specified for the design of 
the ideal (or loss-free) supersonic nozzle except for VE, the Prandtl -Meyer angle cor- 
responding to the exit Mach nximber, ME. This output, also, contains the nozzle sub- 
sonic section coordinates. 

(2) Output 2 gives a table of the coordinates (X and Y) of the supersonic portion of 
the ideal nozzle based on a throat half-height of unity together with the corresponding 
Mach niimber (M) and total to static pressure ratio (FTPS). 
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(3) Output 3 gives the input to the boundary-layer program. The first three lines of 
this output are a listing of the input data except for NST, which is the total number of 
points. The remaining output is the nozzle coordinates for the suction surface which 
contains 10 points on the straight section: 

X X- coordinate, m (ft) 

Y Y-coordinate, m (ft) 

POPTZ static to total pressure ratio corresponding to X and Y 

TWAL wall temperature, K (°R) (wall temperatures are assumed equal to total 
temperature, TTZ) 

(4) Output 4 lists the properties calculated for the upstream conditions: 

2 2 

PSZ upstream static pressure, N/m (Ib/ft ) 

TSZ upstream static temperature, K (°R) 

UZ upstream velocity, m/sec (ft/sec) 

ASZ (ATZ) speed of soimd based on upstream static (total) temperature, m/sec 

(ft/sec) 

RHSZ (RHTZ) static (total) density based on upstream static (total) temperature, 
kg/m^ (slug/ft^) 

MUSZ (MUTZ) dynamic viscosity based on upstream static (total) temperature, (N) 

(sec)/m^ ((lbf)(sec)/ft^) 

NUSZ (NUTZ) kinematic viscosity based on upstream static (total) temperature, 
m^/sec (ft^/see) 

CP specific heat at constant pressure, J/(kg)(K) ((ft)(lbf)/(slug)(°R)) 

PR Prandtl number 

TC thermal conductivity, J/(m)(sec)(K) ((ft)(lbf)/(ft)(sec)(°R)) 

ARCL total distance along surface in x-direction, m (ft) 

The next part of this output gives the variables describing the flow along the surface: 
PRES static pressure, N/m^ (Ibf/ft^) 

UE free-stream velocity, m/sec (ft/sec) 

ME free-stream Mach number 

POPTZ static to total pressure ratio 

VOVCR ratio of the velocity to the critical velocity 

(5) Output 5 corresponds to KMAIN. It indicates the regions of laminar and turbu- 
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lent boundary layers, and the stations at which instability, transition, and separation 
occur. It gives all the principal boimdary- layer output parameters: 

X X- coordinate, m (ft) 

S surface length, m (ft) 

DELSR displacement thickness, m (ft) 

THET momentum thickness, m (ft) 

DELTA boundary -layer thickness, m (ft) 

FORM compressible form factor 

FORMI incompressible form factor 

The next part of this output gives the skin-friction and heat-transfer parameters: 

CF skin-friction coefficient at the wall 

TAUW shear stress at the wall, N/m^ (Ibf/ft^) 

RTH momentum-thickness Reynolds number 

DTDY slope of temperature profile at the wall, K/m (°R/ft) 

NUSS local Nusselt number 

HTRAN heat transfer per unit area, J/(sec)(m^) ((ft)(lbf)/(sec)(ft^)) 

CRN Reynolds analogy parameter 

(6) Output 6 lists the aftermixing properties. Variable names ending in 0 refer to 
the upstream station, those ending in 1 refer to the nozzle exit (before mixing), and 
those ending in 2 refer to the mixed conditions downstream of the nozzle. The output 
consists of nozzle exit free-stream Mach number (XMFSl), nozzle exit spacing 
(SPACING) in meters (ft), traiUng-edge thickness (TE) in meters (ft), Mach numbers 
(XMl and XM2), axial Mach numbers (XMXl and XMX2), critical velocity ratios 
(V/VCRl and V/VCR2), flow angles measured from axial direction (ALPHl and ALPH2) 
in degrees, total to total pressure (PT2/PT0), total to static pressure ratio (PT0/P2), 
static to total temperature ratio (T2/TT0), nozzle kinetic energy loss coefficient 
(EBAR2), and nozzle efficiency (ETA-N). For the case with boundary- layer correction, 
the displacement thickness (DELSR) and the momentum thickness (THET) on the suction 
surface calculated in AFMIX are printed out. The (N+2) values were used in calculating 
the aftermixing properties. These values were calculated because the length of the 
straight section on the suction surface increases when the boundary-layer correction is 
added and a constant nozzle angle is desired. 

(7) The final output 7 lists the coordinates of the supersonic portion of the nozzle 
corrected for the boimdary- layer displacement thickness. The values of the nozzle 
coordinates (X and Y) are with the nozzle axis axial, while XTU, YTU, XTL, and YTL 
are the coordinates of suction and pressure surfaces, respectively, with the axis at the 
nozzle flow angle (see fig. 1). 

10 



TABLE H. - SAMPLE OUTPUT 


TWO DIMENSIONAL SUPERSONIC NOZZLE WITH A SHARP-EDGED THROAT 
"VE 19.999 ME = 1.770 GAMMA => 1.390 DELTA V =0.10000 YTH = O.llOOOE-01 ALPl = 20.000 

SUBSCNIC SECTION (INPUT) 


X Y PSPT 


0. 

5.82180 

0.99303 


10.19990 

4.18240 

0.98640 


14.«727C 

3.46350 

0.98003 


17.78360 

2.96350 

0.97250 


20.C6510 

2.59680 

0.96383 


22.29AA0 

2.23850 

0.95060 


23.55990 

2.03510 

0.93947 


25.42200 

1.73580 

0.91433 


26.32850 

1.59010 

0.89562 


27.20830 

1.44870 

0.87027 


27.88580 

1.33980 

0.84302 


28.96970 

1.16560 

0.77164 


29.57630 

1.06810 

0.69507 


29.87120 

1.02070 

0.62351 


30.00000 

1.00000 

0.52828 


SUPERSCNIC 

SECTION 



X 

Y 

M 

PTPS 

0. 

1.00000 

1.00000 

1.8868924 

1.12958 

1.19752 

1.44655 

3.3859348 

1.34182 

1.23233 

1.46352 

3.4692084 

1.50189 

1.25713 

1.48044 

3.5546543 

1.64083 

1.27741 

1.49731 

3.6423670 

1.76662 

1.29492 

1.51414 

3.7324436 

1.88990 

1.31046 

1.53095 

3.8249831 

2.00713 

1.32444 

1.54773 

3.9200864 

2.12195 

1.33712 

1.56450 

4.0178579 

2.23541 

1.34864 

1.58126 

4.1184044 

2.34826 

1.35911 

1.59801 

4.2218362 

2.46108 

1.36859 

1.61477 

4.3282676 

2.57430 

1.37710 

1.63154 

4.4378166 

2.68829 

1.38467 

1.64833 

4.5506049 

'2.80334 

1.39130 

1.66513 

4.6667590 

2.91972 

1.39699 

1.68196 

4.7864097 

3.03764 

1.40173 

1.69881 

4.9096938 

3.15732 

1.40548 

1.71570 

5.0367521 

3.27894 

1.40824 

1.73262 

5.1677312 

3.40266 

1.40997 

1.74958 

5.3027837 

3.52867 

1.41063 

1.76659 

5.4420699 

3.52867 

1.41063 

1.76659 

5.4420699 



table n. - Continued. SAMPLE OUTPUT 


BOUNDARY LAYER SECTION 


R 

PTZ 

TTZ 

UPMACH 

1716.00 

12696.00 

1010.00 

0. 

NST 

NVP 

NTURB 

KEM 

45 

5 

-0 

-0 

OLAR 

TLAM 

DTURB 

TTUR8 

0.000250 

0.000100 

-0. 

-0. 

KPRE 

KGRAO 

KSDE 

KLAM Kt 

-0 

-0 

-0 

-0 

X 

Y 

POPTZ 

TWAL 

0. 

0.06404 

0.993030 

1010. COOO 

0.11220 

0.04601 

0.986400 

1010.0000 

0.16140 

0.03810 

0.980030 

1010. COOO 

0.19562 

0.03260 

0.972500 

1010. COOO 

0.22072 

0.02856 

0.963830 

1010. COOO 

0.24524 

0.02462 

0.950600 

1010. COOO 

0.25916 

0.02239 

0.939470 

1010.0000 

0.27964 

0.01909 

0.914330 

1010. COOO 

0.28961 

0.01749 

0.895620 

1010. COOO 

0.29929 

0.01594 

0.870270 

1010.0000 

0.30674 

0.01474 

0.843020 

1010.0000 

0.31867 

0.01282 

0.771640 

1010. COOO 

0.32534 

0.01175 

0.695070 

1010.0000 

0.32858 

0.01123 

0.623510 

1010. COOO 

0.33000 

0.01100 

0.529972 

1010. COOO 

0.34243 

0.01317 

0.295339 

1010. COOO 

0.34476 

0.01356 

0.288250 

1010. COOO 

0.34652 

0.01383 

0.281321 

1010.0000 

0.34805 

0.01405 

0.274547 

1010. COOO 

0.34945 

0.01424 

0.267921 

1010.0000 

0.35079 

0.01442 

0.261439 

1010. COOO 

0.35208 

0.01457 

0.255096 

1010. COOO 

0.35334 

0.01471 

0.248889 

1010. COOO 

0.35459 

0.01484 

0.242812 

1010. COOO 

0.35583 

0.01495 

0.236864 

1010. COOO 

0.35707 

0.01505 

0.231039 

1010.0000 


KPROF 

-0 



PPELININARY CALCULATIONS 


PSZ 

= 

1Z696 . 00000 

TSZ 


1010.0000 

UZ 

as 

0. 

ASZ 


1552.1251 

ATZ 

- 

1552.1251 

RHSZ 


0.73253A8E-02 

RHTZ 

s 

0.7325348E-02 

POSZ 


C.6048975E-06 

KUT5Z 


0.6048975E-06 

NUSZ 

3; 

0.8257594E-04 

NUTZ 


0.8257594E-04 

CP 


6115.9999 

PR 


0.66353 

TC 

3 

0.5663117E-02 

ARCL 


0.4586 


STATfIGN 

PRES 

UE 

ME 

POPTZ 

VOVCR 

•1 

12607.50879 

155.63106 

0.100368 

0.993030 

0.109611 

2 

12523.33435 

217.65661 

0.140501 

0.986400 

0.153295 

3 

12442.46082 

264.05568 

0.170607 

0.980030 

0.185974 

4 

12346.85986 

310.29330 

0.200699 

0.972500 

0.218539 

5 

12236.78564 

356.43050 

0.230830 

0.963830 

0.251034 

'6 

12066.81750 

417.57410 

0.270953 

0.950600 

0.294097 

7 

11927.51111 

463.19551 

0.301052 

0.939470 

0.326228 

8 

11608.33362 

553.70406 

0.361250 

0.914330 

0.389973 

9 

1137C.79150 

613.41231 

0.401368 

0.895620 

0.432026 

10 

11046.94788 

687.30137 

0.451530 

0.870270 

0.484066 

11 

10702.98181 

760.24225 

0.501683 

0.843020 

0.535438 

12 

9796.74146 

930.96562 

0.622016 

0.771640 

0.655678 

13 

8824.60864 

1094.82568 

0.742300 

0.695070 

0.771085 

14 

7916.08295 

1238.37828 

0.852526 

0.623510 

0.872189 

15 

6728.52344 

1419.85098 

1.000000 

0.529972 

1.000000 

16 

3749.62912 

1892.13937 

1.446551 

0.295339 

1.332632 

17 

3659.62442 

1907.82626 

1.463524 

0.288250 

1.343681 



TABLE n. - Continued. SAMPLE OUTPUT 

PRIIICIP4I. BOUNDARY LAYER INFORMATION 

INSTABILITY DOES NOT OCCUR 
TRANSITION DOES NOT OCCUR 
SEPARATION DOES NOT OCCUR 

LAMINAR BOUNDARY LAYER - STATIONS 1 TO 45 
TURBULENT BOUNDARY LAYER DOES NOT OCCUR 


STATION 

X 

S 

DELSR 

THET 

DELTA 

FORM 

FORMI 

I 

0. 

0. 

0.000243 

0.000100 

0.000832 

2.4340 

2.4276 

2 

0.112199 

0.113639 

0.00025C 

0.000100 

0.000921 

2.5054 

2.4927 

3 

0.161400 

0.163471 

0.000232 

0.C00092 

0.000886 

2.5356 

2.5167 

4 

0.195620 

0.198130 

0.000209 

0.000082 

0.000815 

2.5593 

2.5330 

5 

0.220716 

0.223549 

0.000188 

0.000073 

0.000740 

2.5781 

2.5432 

S 

0.245238 

0.248386 

0.000166 

0.000064 

0.000656 

2.6032 

2.5549 

7 

0.259159 

0.262485 

0.000152 

0.000058 

0.000602 

2.6209 

2.5613 

8 

0.279642 

0.283231 

0.000131 

0.000049 

0.000520 

2.6641 

2.5778 

9 

0.289613 

0.293330 

0.000120 

0.000045 

0.000474 

2.6907 

2.5839 

10 

0.299291 

0.303132 

0.000110 

O.CC0040 

0.000441 

2.7452 

2.6090 

11 

0.306744 

0.310681 

0.000100 

0.000036 

0.000385 

2.7623 

2.5950 

12 

0.318667 

0.322756 

0.000089 

0.000030 

0.000380 

2.9886 

2.7218 

13 

0.325339 

0.329515 

0.000072 

0.000024 

0.000244 

2.9676 

2.6006 

14 

0.328583 

0.332800 

0.000077 

0.000021 

0.000412 

3.6684 

3.1119 

15 

0.330000 

0.334235 

0.000068 

0.000018 

0.000314 

3.7926 

3.0408 

lf6 

0.342425 

0.346849 

0.000089 

0.000022 

0.000210 

4.0912 

2.6696 

17 

0.344760 

0.349215 

0.000096 

0.000023 

0.000230 

4.1602 

2.6946 

18 

0.346521 

0.350997 

0.000102 

0.000024 

0.000245 

4.2218 

2.7138 

19 

0.348049 

0.352541 

0.000107 

0.000025 

0.000257 

4.2789 

2.7295 

20 

0.349455 

0.353960 

0.000111 

0.000026 

0.000268 

4.3328 

2.7426 

31 

0.350789 

0.355305 

0.000116 

0.000026 

0.000279 

4.3848 

2.7539 

22 

0.352078 

0.356604 

0.000120 

0.000027 

0.000289 

4.4353 

2.7638 

23 

0.353341 

0.357874 

0.000124 

0.000028 

0.000298 

4.4848 

2.7727 

ZA 

0.354590 

0.359129 

0.000128 

0.000028 

0.000307 

4i5337 

2.7807 

25 

0.355831 

0.360376 

0.000132 

C. 000029 

0.000316 

4.5822 

2.7882 

Zb 

0.357072 

0.361621 

0.000135 

0.000029 

0.000324 

4.6304 

2.7951 

27 

0.358317 

0.362870 

0.000139 

0.000030 

0.000333 

4.6786 

2.8015 



ST«TION 

CF 

TAUW 

RTH 

DTOY 

NUSS 

HTRAN 

CRN 

1 

0. 

0. 

187.2 

0. 

0. 

0. 

2.799 

2 

0.0C287 

0.49111 

259.9 

-978.37 

154.81 

-5.5407 

4.842 

3 

0.00277 

0.69284 

287.2 

-1340.26 

207.28 

-7.5900 

6.050 

4 

0.0C276 

0.94582 

298.9 

-1911.08 

259.42 

-10.8227 

6.806 

5 

0.0C276 

1.23823 

304.0 

-2708.88- 

314.43 

-15.3407 

7-221 

6 

0.0C279 

1.69433 

305.9 

-4098.73 

385.15 

-23.2116 

7.650 

7 

0.0C281 

2.07543 

305.2 

-5445.08 

439.44 

-30.8361 

7.824 

8 

0.0C290 

2.97685 

300.7 

-8772.62 

534.61 

-49.6803 

8.327 

9 

0.0C293 

3.61860 

296.3 

-11754.33 

604.46 

-66.5661 

8.369 

10 

0.0C311 

4.67552 

289.5 

-15313.05 

648.22 

-86.7196 

9.300 

11 

O.OC300 

5.35578 

281.7 

-22139.43 

785.05 

-125.3782 

8.151 

12 

0.0C399 

9.77562 

259.1 

-29728.16 

730.29 

-168.3540 

13.567 

13 

0.0C322 

9.81642 

224.3 

-78886.76 

1430.57 

-446.7449 

6.038 

14 

0.0C809 

28.34312 

197.3 

-50701.23 

725.80 

-287.1270 

30.681 

15 

0.0C817 

31.98068 

162.6 

-80457.29 

879.94 

-455.6390 

25.012 

16 

0.0C259 

10.03554 

146.8 

-332898.43 

2127.49 

-1885.2427 

2.528 

17 

0.CC268 

10.28577 

154.4 

-298177.55 

1887.17 

-1688.6143 

2.917 

18 

0.0C274 

10.45993 

158.5 

-274410.61 

1717.63 

-1554.0193 

3.249 

19 

0.00280 

10.57238 

161.2 

-256461.62 

1587.04 

-1452.3721 

3.542 

20 

0.00284 

10.63108 

163.0 

-242547.11 

1-483.79 

-1373.5726 

3.795 

21 

0.CC287 

10.65249 

164.3 

-231373.06 

1399.37 

-1310.2926 

4.017 

22 

0.0C290 

10.64292 

165.3 

-222221.00 

1328.97 

-1258.4635 

4-210 

23 

0.0C291 

10.60644 

166.0 

-214596.34 

1269.22 

-1215.2841 

4.376 

24 

0.0C293 

10.54865 

166.5 

-208095.97 

1217.46 

-1178.4717 

4.521 

25 

0.0C294 

10.47273 

166.8 

-202460.90 

1171.94 

-1146.5597 

4.646 

26 

O.OC295 

10.38160 

167.1 

-197493.29 

1131.32 

-1118.4275 

4.754 

27 

0.0C295 

10.27764 

167.2 

-193047.90 

1094.63 

-1093.2528 

4.848 

28 

C.0C295 

10.16284 

167.2 

-189014.53 

1061.12 

-1070.4113 

4.929 

29 

0.0C295 

10.03898 

167.1 

-185306.65 

1030.23 

-1049.4132 

4.999 

30 

0.0C295 

9.90755 

167.0 

-181858.09 

1001.49 

-1029.8836 

5.060 

31 

0.0C294 

9.76977 

166.9 

-178617.79 

974.57 

-1011.5334 

5.113 

32 

0.0C291 

9.54961 

166.7 

-172338.05 

931.84 

-975.9705 

5.213 

33 

0.0C291 

9.40499 

166.5 

-168505.05 

903.13 

-954.2638 

5.283 

34 

0.0C286 

9.13631 

166.2 

-173826.39 

923.69 

-984.3992 

5.006 

35 

0.0C170 

5.35489 

165.9 

-250520.49 

1320.16 

-1418.7268 

2.048 

36 

0.0C144 

4.53268 

184.1 

-227632.28 

1226.94 

-1289.1082 

1.908 

37 

0.0C138 

4.34516 

202.1 

-206044.66 

1135.38 

-1166.8550 

2.021 

38 

0.0C128 

4.01812 

218.5 

-190536.64 

1072.85 

-1079.0312 

2.021 

39 

0.00119 

3.75531 

233.8 

-178074.28 

1024.11 

-1008.4554 

2.021 

40 

0.0C112 

3.53813 

248.2 

-167775.60 

985.08 

-950.1328 

2.021 

41 

0.0C107 

3.35474 

261.7 

-159079.22 

953.16 

-900.8842 

2.021 

42 

0.0C102 

3.19719 

274.6 

-151608.55 

926.64 

-858.5769 

2.021 

43 

0.OC097 

3.05994 

286.9 

-145100.33 

904.33 

-821.7201 

2.021 

44 

0.0C093 

2.93898 

298.8 

-139364.19 

885.35 

-789.2356 

2.021 

45 

0.0C090 

2.83131 

310.1 

-134258.63 

869.07 

-760.3223 

2.021 



TABLE n. - Concluded. SAMPLE OUTPUT 


XMFSl 

ALPHl 

XPFSl 

ALPHL 


AFTERMIXING PROPERTIES 
NOZZLE WITH NO BOUNDARY LAYER CORRECTION 


=1.7666 

SPACING 

=.090740 TE =0. 

XM2 =1.7532 

V/VCRl 

= 1.523 XMXl 

= 0.604 

XMX2 = 

0.593 


i 70.001 

ALPH2= 

70.241 PT2/PT0= 0.9861 

PT0/P2= 

5.407 

T2/TT0= 0.6253 

V/VCR2= 

1.515 

EBAR2=0. 00651 

ETA-\I=0.9935 

NOZZLE WITH 

BOUNDARY LAYER CORRECTION 







=1.7666 

SPACING 

=.092147 TE =0. 

XH2 =1.7533 

V/VCRl 

= 1.523 XMXl 

= 0.604 

XMX2 = 

0.593 


69.999 

ALPH2= 

70.237 PT2/PT0= 0.9862 

PT0/P2= 

5.407 

T2/TT0= 0.6252 

V/VCR2= 

1.516 

EBAR2=0. 00643 

ETA-N=0.9936 


DELSRIN-i-l) S0.C00315 DELSRIN + 2) =0.000315 THET(N+1) =0.000065 THETIN+Z) =0.000065 



r 


NOZZLE COORDINATES SUPERSONIC SECTION 


NCZZLE WITH BOUNDARY LAYER CORRECTION 


X 

Y 

XTU 

Q. 

0.011068 

-0.010400 

0*012425 

0.013261 

-0.008211 

C. 014760 

0.013652 

-0.007780 

0.016521 

0.013931 

-0.007439 

C. 018045 

0.014159 

-0.007131 

0.019455 

0.014356 

-0.006835 

0.020789 

0.014531 

-0.006543 

0.022078 

0.014689 

-0.006251 

0.023341 

0.014832 

-0.005953 

0.024590 

0.014963 

-0.005649 

0.025831 

0.015082 

-0.005337 

0.027072 

0.015190 

-0.005014 

0.028317 

0.015288 

-0.004680 

0.029571 

0.015375 

-0.004333 

0.030837 

0.015452 

-0.003972 

0.032117 

0.015519 

-0.003597 

0.033414 

0.015575 

-0.003207 

0.034731 

0.015621 

-0.002799 

0.036068 

0.015656 

-0.002375 

0.037429 

0.015679 

-0.001931 

0.038815 

0.015685 

-0.001462 

0.047342 

0.015703 

0.001438 

0.055869 

0.015721 

0.004337 

0.064396 

0.015738 

0.007238 

0.072922 

0.015753 

0.010140 

0.081449 

0.015768 

0.013043 

0.089976 

0.015782 

0.015946 

C.0985C3 

0.015795 

0.018851 

0.107030 

0.015807 

0.021756 

0.115556 

0.015819 

0.024661 

0.124083 

0.015830 

0.027567 

0.125000 

0.015832 

0.027879 

0.125405 

0.015832 

0.028017 


YTU 

XTL 

YTL 

0.003786 

0.010400 

-0.003786 

0.016212 

0.016712 

0.007140 

0.018539 

0.017877 

0.009200 

0.020289 

0.018741 

0.010759 

0.021803 

0.019478 

0.012118 

0.023192 

0.02C144 

0.013371 

0.024505 

0.020765 

0.014565 

0.025771 

0.021354 

0.015723 

0.027007 

0.021921 

0.016860 

0.028224 

0.022471 

0.017988 

0.029432 

0.023007 

0.019114 

0.030635 

0.023534 

0.020243 

0.031838 

0.024051 

0.021380 

0.033046 

0.024562 

0.022529 

0.034262 

0.025068 

0.023691 

0.035488 

0.025568 

0.024871 

0.036726 

0.026065 

0.026071 

0.037979 

0.026558 

0.027292 

0.039248 

0.027049 

0.028538 

0.040535 

0.027536 

0.029809 

0.041839 

0.028015 

0.031109 

0.049858 

-O.OOCOOO 

-0.000000 

0.057876 

-0.000000 

-0.000000 

0.065895 

-0.000000 

-0.000000 

0.073912 

-0.000000 

-0.000000 

0.081930 

-0.000000 

-0.000000 

0.089947 

-0.000000 

-0.000000 

0.097964 

-0.000000 

-0.000000 

0.105981 

-0.000000 

-0.000000 

0.113997 

-O.OOCOOO 

-0.000000 

0.122013 

-O.OOCOOO 

-0.000000 

0.122875 

-O.OOCOOO 

-0.000000 

0.123256 

-0.000000 

-0.000000 



PROGRAM DESCRIPTION 


The program SSN is divided into two main parts, subroutines NOZZL and BL. This 
division was necessary because of the storage requirements of the IBM 7094/7044 direct 
coupled system which has a 32 768 word core (77777g). An overlay system was used 
with SSN as the common link between NOZZL and BL. NOZZL requires 51216^g^ words 
of storage and BL requires 43535^g^ words of storage. SSN uses 14134^g^ words of 
storage. All quantities calculated in NOZZL and needed by BL are stored in SSN 
through two common blocks and the CALL statement. 


Subroutine NOZZL 

This subroutine designs a loss-free supersonic nozzle by the method of character- 
istics. The description of the method, program, and main dictionary of variables is 
given in reference 1. The length of the straight section on the suction surface of the 
nozzle is determined by the nozzle angle. 

The calling sequence for NOZZL is as follows: 

CALLNOZZL(X, Y, POPTZ, SP; NST, NCP, ALPl, GAM, NSP) 

where 

X X- coordinate of nozzle contour point (output) 

Y Y-coordinate of nozzle contour point (output) 

POPTZ static to total pressure ratio (output) 

SP spacing (output) 

NST total number of points (output) 

NCP number of points excluding the ten points for the straight section (output) 

ALPl nozzle flow angle at station 1, measured from tangential direction (output) 

GAM specific heat ratio (output) 

NSP number of subsonic points including throat point (output) 
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The program variables for NOZZL that are not in reference 1 are 


C 

CH 

M 

PSPT (NC) 
PTPS (NC) 
S 
Z 


chord based on throat haH- height of unity 
chord, m (ft) 

Mach number 

static to total pressure ratio at a contour point 
total to static pressure ratio at a contour point 
spacing, m (ft) 

length of straight section, m (ft) 


Subroutine BL 

This subroutine calculates the bomdary- layer characteristics, mixing losses, and 
corrected nozzle contoxir. The major subroutines for BL are INPUT, PRECAL, 
LAMNAR, TURBLN, PROFIT, AFMIX, and NOZZLC. These, in turn, call several 
other subroutines. All these subroutines and their relation are shown in figure 3. 



Figure 3. - Calling relation of BL subroutines. 


Subroutines INPUT, PRECAL, LAMNAR, TURBLN, and PROFIT . - These subrou- 
tines and their auxiliary subroutines shown in figure 3 calculate the boundary- layer 
characteristics. The methods used for calculating the boundary- layer characteristics 
and the program description are given in reference 3. 

Since a highly favorable presstme gradient exists, the range of the equations of ref- 
erence 5 was extended by the method given in reference 5. These changes were made in 
LAMNAR and are for SW (temperature function at the wall) = 0. The changes are noted 
in the program with comment cards. Two curve fit ranges were also extended. RCRIT 
and DIFF were extended as follows: 
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RCRIT = 8. 3163, when SHARK is greater than 0. 07 

DIFF = 44 000 KBAR + 700, when KBAR is greater than 0. 03 

Subroutine AFMIX . - Subroutine AFMIX takes the boundary- layer parameters and 
free-stream conditions at the nozzle exit and calculates the aftermixing conditions by 
the method given in reference 4. In this loss model, the flow sufficiently downstream of 
the blade row is assiimed to be mixed to a uniform condition. The aftermixing condi- 
tions are calculated for two cases - nozzle without boundary- layer correction, and noz- 
zle with boundary -layer correction. 

For the case with no boimdary- layer correction, the displacement and momentxun 
thicknesses on the suction and pressure surface and the free-stream conditions at noz- 
zle exit (station 1) are used to calculate the aftermixing conditions (station 2). 



Figure 4. - Quantities required for suction surface calculations at nozzle exit. 


For the case with boundary -layer correction, the displacement and momentum 
thicknesses on the suction surface at the nozzle exit must be calculated. When the dis- 
placement thickness is added to the nozzle profile, the length of straight section on the 
suction surface must be increased to keep the nozzle angle constant. Therefore, from 
figure 4, 


( 1 ) 

The change in length of straight section is 

AZ = Z - (X - X ) (2) 

' n np'' ' ' 

The coordinates of the straight section at this point are 
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(3a) 


X 


n+1 




AZ 


'n+1 



(3b) 


The increase in momentum thickness and displacement thickness is approximately linear 
with distance on the straight section. Therefore, the displacement angle is 


FA =: arctan 



The displacement thickness is 


= 6^ + AZ tan FA 


Similarly, 


FA0 = arctan 



(4) 


(5) 


(6) 


Therefore, the momentum thickness is 

^n+1 = 

As can be seen from figure 4, the length of the straight section has to be increased fur- 
ther to account for the boundary- layer addition. The angle at the nozzle exit at suction 
surface is 


- 180° - 0-1 - (90° + FA) 

W X 


( 8 ) 


From the law of sines 


As _ ^n+1 ~ ^np 
sin (90 + FA) sin 

w 


The change in spacing is 
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As = (8 , - 8 ) sin (8° + ^A) 

"+1 "P sin » 


w 


( 9 ) 


The corrected spacing is 


Sg = s + As 


(10) 


The change in axial length is 


AZj = 


sin O'- 


sin 0 !, 


K - «np) 


w 


( 11 ) 


Therefore, the displacement and momentum thicknesses on the suction surface at the 
nozzle exit are 


^2 = Vl + ^^f (12) 

- ®n+l Is-n FA0 (13) 

The coordinates of the straight section at this point are 


V2 = Xn+1 - ^Zf 


V - V 
n+2 np 


(14a) 

(14b) 


These values of the displacement and momentum thicknesses are used for the suction 
surface in the calculation of the mixing losses. 

The calling sequence for AFMIX is 


CALLAFMIX (ALPHl, TE, SP, ME(NST), NP) 


where 


ALPHl nozzle flow angle at station 1 , measured from axial direction (input) 
TE trailing-edge thickness (input) 

SP spacing (input) 

ME(NST) free-stream Mach number at nozzle exit, station 1 


22 



NP number of points to end of bell portion of the nozzle (last point on pressiire 

surface) 

The program variables for AFMIX are 

A 1 - 6* - - 0* 

A1 1 - 6* - 6^^ 

AFSl eq. (B4) of ref. 4 

ALPHl see input 

ALPW Q!,„, fig. 4 

w 

C eq. (C16) of ref. 4 

D eq. (C18) of ref. 4 

DELP displacement thickness on pressure surface 

DELS displacement thickness on suction surface 

DELSR displacement thickness 

DELSRT a*, eq. (18a) of ref. 4 

DELZ AZ 


DELZF 

DENCH 

DENSRl 

DENSR2 

DS 

DF 

DTE 

EBAR 

ETAN 

FA 

FAMT 

GAM 

KODE 

N 


AZj 

part of eq. (C21) of ref. 4 
(P/P’)j 

eq. (C21) of ref. 4 
change in spacing, As 
change in surface length 
a^g, eq. (1) of ref. 4 

nozzle kinetic energy loss coefficient, eq. (24) of ref. 4 
nozzle efficiency 

displacement thickness angle, eq. (4) 
momentum thickness angle, eq. (6) 
specific heat ratio 

program control for different solutions 
NST 


NP see input 
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PR2 (p/p*)2 

PT2PTO eq. (C22) of ref. 4 

PTOP2 p;/p2 

SP spacing, s 

THET momentum thickness 

THETA 0* , eq. (18b) of ref. 4 

THETS momentum thickness on suction surface 

THETP momentum thickness on pressure surface 

VXVCR2 ratio of axial velocity to critical velocity at station 2, eq. (C20) of ref. 4 
XX S cos 0!j 

XXX X- coordinate 

YYY Y- coordinate 

Z Z, length of straight section 

Subroutine NQZZLC . - This subroutine calculates the corrected supersonic nozzle 
coordinates by adding the displacement thickness to the ideal nozzle coordinates. The 
coordinates are calculated for the nozzle with the nozzle axis horizontal. The coordi- 
nates are also calculated for the nozzle with the nozzle axis at the nozzle angle. 

The calling sequence for NOZZLC is 

CALL NOZZLC (NSP, ALPHl, NP) 

where 

NSP munber of subsonic points including throat point 

ALPHl nozzle angle 

NP last point on pressure surface 

The program variables are 

XSUB length of subsonic portion, m (ft) 

XXT X- coordinate of the suction surface (nozzle axis at nozzle angle), m (ft) 

XXTL X-coordinate of pressure surface (axis at angle), m (ft) 

XXX X-coordinate (nozzle axis horizontal), m (ft) 

YYC Y- coordinate of suction surface (axis horizontal), m (ft) 

YYCL Y- coordinate of pressure surface (axis horizontal), m (ft) 
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JJ f' 


YYT Y-coordinate suction surface (axis at angle), m (ft) 
YYTL Y-coordinate pressure surface (axis at angle), m (ft) 


PROGRAM LISTING 


ilR^TC SSN OfcCK 

C StJPFRSnNIC NU77LE 

COMMON /Cl/ GAM. K .PT/ .TTZ . UPMACW ,MST ,NVP , NTUKB, KPVM ,KEM, KS MT H , 
IKSPLN.M E.KATCH.CTHFT.DLAM.TLAM.lJTURa.TTORB.KPRE.KGRAD.KSDE .KLAM, 
7KMAIN.KPRUE,X< 100) ,Y( 100 ) .PRES (100 ) .UE (100) . ME ( 100 ) , POPT Z ( 1 DO ) . 
3VOVf,R( 100 ) . TRAL( 100) 

COMMON/CP/P S /, T SZ. UZ. ASZ .AIZ.RHSZ.RHT Z.MUSZ. MUTZ.NUSZ. NjrZ, CP, 
IPR.TC.ARCL 
REAL ME 

1 CALL i\l[)ZZL< X.Y.PUSTZ.SP .NST.N: P.ALPl.&AM.NSP) 

CALL RL ( SP .NCP .ALP 1 .MSP) 

GO TO 1 
END 


SIRFTC NOZZLE UECK 

•SUBRnuriNE NOZZUXCP.YCP.PSPT.S .NOP.NCP.ALPl .GAM.NS) 

C TrO dimensional supersonic nuzzle with a SHARP-EDGED THROAT 

DIMENSION XC(1500) .YC (1500 .PHK1500) .U(3DDD) .MODDO). 

IX ( 1500.2) , Y( 1500.2) . PIPS (3000) .XC P ( 100 ) , YCP( 1 00) .PSPT(IDD) 
real m. mc . me .mec. mt 
REAO( 5.3) ME .0 V .GAM. NP . N S. ALPl ,YTH 
FORMAT ( 3F6.3.2I 2,EC.3,E10. 5) 

READ( 5,9) ( XCP( I ) ,YCP( I ) ,PSPT( I ) ,1=1 ,NS) 

9 F0RMAT(3F10.1) 

NCP=NS+1 

ALPl=ALPl=f‘-017453 

VE = SUrT( (GAM+1.0) / (GAM-1. 0) ) * A TAN (SORT ( (GAM-l. 0 ) /'( G AM+ 1. 3 ) 
l*(MF«=*2-1.0) )) - ATAN(S0HT(ME**2-l.(j) ) 

0V= UY# .017 453 
XT=0.0 
YT=1.0 
MT =1.0 

PTPST = ( (GAM+1.0)/2.0)**(GAM/(GAM-1.0)) 

XMAX = INT( . 5*VE/3V + 1.5) 

0V=VF/(2.O*FLOAT(KMAX-l) ) 

NMAX = KM AX 
WKlTE(tt,251 ) 
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251 FORMATt IH l,3tX, 59HTisD niMtNSlQNAL SUPERSONIC NOZ2LE WITH A SHARP-E 
lOfiEO THROAT} 

ve= WE* 57 -29 58 
OV =OW*57,29 58 
ALPl = ALP 1*57. 2958 

WRITE<6,252) WE. ME , GAM .0 W , YTH, ALP I 

252 FORMAT! 1H0.4HWE =. F8. 3 . 10X.4HME = ,F 8. 3 , 10 X , 7HGAMMA =,F6.3,iaX, 
19HDELTA W =,F7.5. 8X.5HYTH =. t 12. 5 . 8X .6HALP 1 =,F7.3} 

WR ITE(6.247 ) 

247 FOkMATl IHU, 5X, 24HSUBSUNIC SECTION (INPUT)) 

MRITE(6.248 ) 

248 FORMAT! 1H0.5X, IHX, 14X.1HY . 14X.4HPSPT) 

WRITE! 6.2 55) ! XCP ! I ) . YC P ! I ) .P SPT ! I ) .1 =l .NS) 

255 FORMAT! IH . F 10. 5. 2F 15. 5) 

WR ITF!6.246) 

246 FORMAT! IHO. 5 X. 18HSUPERS0NIC SECTION) 

WRITF!6.253> 

253 FORMAT! 1H0,5X, IHX. 14X , IH Y , 14X . IHM. 14X , 4HP TPS ) 

WR ITE! 6.254)XT« YT.MT. PTPST 

A!P1 = ALPl* .017453 
0W= DW* .0 174 53 
C ROW ANGLE CALCULATION 

DO 5 K=1.KMAX 

5 PHIIK) = FLUATlK-l )*UW 

C MACH ANGLE CALCULATION 

NK = 2*KMaX-1 

00 6 1=1, NK 

6 CALL UA! I ,UW .GAM.U! 1 ) ) 

C REGION NEAR TEE THROAT 

DO ICO K=1.KMAX 

NC = I 

N=1 

1 = K«-2*< NC- 1) 

IF IK.NE.I) GO TO 10 

SLOPEl = -TAN! !U! I )*U! 1*1) )/2.0-(PHl (K)+PHI(K<-1) )/2.0) 

X !K,N)= - l.C/SLOPEl 
Y!K,N ) =0.0 
GO TO 100 

10 IF (K.FU.KMAX) GO TO 20 

SLOPEl = -lAM ( U! I ) +U! I + l) )/2.0-!PHI !K) *PHI!K*1) ) /2.0) 

SL0PE2. = TAN! (U! I )+U( 1+ 1) ) /2.U+! PHI (K) +PHI (K-1) )/2.0) 

X!K,N) = (1.0-! Y!K-1,N)-SLUPE2*X!K-1,N) ) ) / 

1! SL0PF2-SLUPF1) 

Y!K.N)=Y! K-l.N) + SLaPE2*(X!K,N)-X(K-l ,N) ) 

(iO TO 100 

20 SL0PE1= TANIPHKK)) 

SL0PF2 = TAN! !U( l)+U( I+i) )/2.C+!PHI (K)*PHI (K-1) )/2.0) 

X(K.N) = ( 1.G-! Y! K-l,N)-SL0PE2*X! K-l ,N) } ) / 

1! SLUPE2-SL0PF1) 

Y ( K,N)=Y! K-l .N ) +SL0PE2*1X! K.N) -X! A-l ,Nl i 
XC! NC )=X! KMAX.N) 

YC(NC)= Y(KMAX.N) 

M! NC) = 1 ,0/SIN!U! I ) ) 

P TPS (NO =! 1.0* ! GAM-1. 0)*M!NC) **2/2.0)**! GAM/ (GAM-1.0) ) 

100 CONTINUE 

IF (MUDINC.NP) .EO.O) WRITE(6,254) XC ( NC) . YC ( NC) , M ! NC) , PT PS ( NC ) 

IF (MUDINC.NP). NE.O) GO TO 203 
XCP(NCP)=XC!NC) 


26 



YCP< NCP l=YC( Nt: 1 
PSPTJNCP) = 1.0/PTPS( NO 
iNi(.P=NCP +1 

C kEUIUN DUiriNSTkEAM UE THE THROAT 
20^ CtlNTJNUF 
20*1 N = 2 

KMAX = Kf'.AX-l 
nC=NL+1 

i)U 200 K = 1,KMAX 
( = K *-2*( NO- IJ 

IF <K.,nF. 1) 0(1 TO 201 

SLOPE! = -TANUUn )+U(l + in/2.O-(PHI(K)+PHX(K+in/2.0) 

X<K ,N »=-( Y< K*1,N-U-SL0PE1*X< K+itN-U l/SLDPEl 
Y ( K ,N ) =0.0 
(ill TO 200 

201 I F< K . EU .<NAX » GO TO 202 

SLOPFl = -TAN(aJUH-U(I + in/2.O-(PHUKH-Pr)I(K+l))/2.0> 

SI OPE? = TANK on ) + U( I*in/2.C+(PHI (K)+PHI(K-1) )/2.0) 

X(K ,N ) = (< Y< K+l.N-U-SLOPEl*X(K+l.N-l) ) - ( Y ( K-1 , N) -S L0PE2 
l*X!K-l.Nn J/(SLnPE2 - SLOPED 
Y( <,NJ=Y< K-1 ,N 1 +SLQPE2* ( X { K . N I -X( K-1 , N) ) 

GO TO 200 

202 SLnPEl=TAN(PHI ( KD 

SL3PE2 = TANUUn l+U(DlD/2.O+4PHUKH-PrU(K-in/2.01 
X<K ,N » = (( Y( K+l.N-D -SLOPE D=X( K + 1,N-1D- (Y( K-1, N)-SL0P£2 
l#X(K-l.Nn 1/ (SL0PE2 - SLOPED 
Y( k.N » = Yt K-1 ,N)+SLnPE2*(X(k.NJ-X4K-l.N ) ) 

XU NC»=X( K.MAX,N I 
YC<NC»=Y( XMAX.NJ 
MINC ) = 1 .0/SIN( U( 1 D 

PTPSl NO » =( 1,0+(GAM-1.0)*M(NC JX=»2/2.0» ♦♦( G AM/ ( GAM-1 .0 ) > 

200 C3MTIMUE 

IF (M(iU(NC,NP» .EU.OJ WRITE(6,254l XC ( NO . YC { NO , M( NO , PTPS ( NO 
IE ( M0U( NC.NP) .NE.O GO TO 249 
XCP( NCP »=XC(NC) 

YLPINCP » = YC(NC J 
PSPT(NGP) = 1.0/PTPS(NO 
N CP =N CP +i 

249 IE (KMAX. EU.2I GO 10 250 
1)0 204 K=1,KMAX 

X(K. 1 ) = X( K. 2 ) 

204 Y<K.D = Y1K.2) 

GO TO 205 

250 NMA=NMAX-1 

lElXCPi NCP-1 ».E0.XC<NMAD GO TO 300 
XCP1NCP>=XC( NMA» 

YCPINCP ) = YC(NMA » 

PSPDnCP » = 1.0/PTPSINMA) 

GO TO 301 
30U NCP=NCP-l 

301 »<kITE (6, 254) XU M«A) ,YC (NMA ) .M(NMA) . PTPS(NMA) 

254 FOkMAKlH . F 10 . 5, F 15 . 5,F 1 5. 5,F 15. 7 ) 

/. =2 .0« YC ( N MA )/ T AN ( AL P 1 ) *Y TH 

S = 2.U<'YC(NMA ) /SIN( Ai.PD #YTH 

C = XCP(nS) +XC(NMA»*-//YTH 

CH = (.*YTH 

PSPT(NS) = 1,0/PTPST 

Nap=NCP+10 

nil 3(. 4 Dl.NOP 
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IF < (r4S+l) ) r,(j TO 302 

XCP( I J = XC.P( I )/U*CH 
YCP( I»=YOP( I J/C*CH 
GO TO 304 

302 IF( I,UE .( NCP + in Gu TO 303 
XtPI n = XCP( I l/C*CH + XCP(NS» 
YCP( I1 = YC.P< I )/C*CH 

GO TO 304 

303 PSPI(I»= PSPKNCP) 

xcpcn =z/io.o*xcpi i-i) 

YCP( I )=YLP( NCP) 

304 CONTIiMOE 
RFTiJRN 
END 


ilbFTO U DECK 

SIJHKOOTINE U/^( I .DV.GAM.U) 

REAL MC.M.MCl 

X2M?S0RT( (GAM+1.0> /(GAM-l.on-.Ol 

V = FI GAT ( I- IJ-i'DV 

IF(V.FO.O.OJ GO TG 5 

EXTERNAL YA 

IF |M(.,NE.1.C) GO TG 7 

DVDMC=.Q2 

GO TG 4 

7 A= SURT( ( GAM*1.0)/(GAM-1,0I ) 
rt= ( GAM-1 .0 )*MC**2-GAM 

C= ( GAM+1 .0 ) /MG*=«=2-GAM 

DyOM(,= Aft( GAM-1. Ot #«C / ( 1 . 0-8**2 )** ( . 5 ) - MG *♦ ( -3 ) ♦ 1 G AM+l . 0 ) / ( l.O- 
1C**2 5 1 

4 MC= M /OVDMC+MC 

IF <MC.lt .X2M1 GO TO 9 
MC=X2M 

WRITE (6.2) MC.NC.K 

2 FORMAT! IHO, 4 IHL IMI I HAS BEEN REACHED. MC SET = TO X2M,5X,4HMC 

1F8.3.5X.4HNC =,I4.5X.3HK =,I4) 

9 Vl= VAI MC ) 

IF IAHS(V-V1 ).LT. .OOCOU GO TO B 
GO TO 7 

8 M = SORT! < <2.0/(GAM4-1.(;))*MC**=2)/(1.0-(<GAM-1.0)/ 

l( GAM+l.O) ♦MC**2) ) ) 

GO TO 6 
<5 M = 1.0 

MC=1.0 

6 0=ARS IN( 1 .0/M) 

k FTURN 
END 
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SIBFTC VAN 


UECK 


FUNCTION VA(XI 

CnMMUN/Cl/UAM,« ,PTZ . TTZ . UPMAC ^ , NST, N VP » NTURB , KPVM ,KEM, K.S MI H, 
IKSPLN.KCr .ICATC^ ,CT HET »DLAM» TLA M» 0 TURB »TT J RB» KPKE t KGRADi <S3E fKLAM» 
2<“|AIN.<PT0F.X( 100 I » YI 100 » .PRE S (I 00) ♦U£(133) » ME (1)3) tPQPTtilD) )* 
^iVOVCkl 100), TWAL(IUC) 

A= SORT((GAM+1.0) /(UAM-1.0) )-1.0 
B= AR.SIN( (GAM-1. 0H‘X**2-GAM) 

C=ARSIN(( (iAM4l,0)/X**2-GAM) 

VA = (3.1415926)/A.0 * A + .5* ( ( A+ 1.0 M'B+C) 

RETURN) 

ENO 


$1E?TC SLAYER DECK 

SUBRCUTINE BL ( SP ,NP t ALP 1 ,NSP ) 

CCMMCN/Cl/GAM,Rt.PTZf TTZ»UPMACH,NSTtNVP,NTURB»KPVMtKEM,KSMTH, 
IKSPLN *KLE , KATCHt CTHET,OLAM, TLAM ,DTUR8 ♦ TTURB.KPRE ,KGRADf KSOE,KLAM, 
2KPAIN»KPR0Ff XdOO ) ,Y ( 100 )» PRES ( 100 ) »UE ( 100 ) tME ( 100 ) ,P0PTZ(100) • 
SVOVCRdOC) ,TWAU 100) 

CCKP0N/C2/PSZ,.TSZ»UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP, 

IPRtlCtARCL 

CCMMCN/C3/X0M(100),Y0M( 100),S( 100 ) , SOL ( 100 ) , AE ( 100 ) »TSE ( 100 ) , 
ITAWLI 100) ,TAWT(100),TBARa00),RW(100)#SW(lC0) »SUTHL(100) t 
2RHSVH 100 ) »RHSE( 100 ) ,HEADW( 100 ) »HEADE( 100) wNUW( 100) »MUBAR( 100) # 

3AA { ICO ) ,eB( 100) »^FF( ICO) » DUOS (100 ),DMOS( 100) ,OMDL( 100) 
COKIfCN/CA/THET( 100) fOELSR( 100) ,DEl,TA( 100) , FORM) 100) t 
IFQRMI dOC ) f FORM! R( 100) »RTH( 100) tRTHI ( 100) tCF ( 100 ) » 
ITAUW1100)»NUSS1100) tOTDY(lOO),HTRAN(lOO)fCRN(lOO) 

CCHM0N/C9/ERR0R* TRANS, SEPRN 
REAL ME 

LOGICAL ERROR, TRANS, SEPRN 
10 CALL INPUT 

REAC 65,99) TE 
99 F0RMAT(F6.1) 

IF (ERROR) GO TO 10 
CALL PRECAL 
IF (ERROR) GO TO 10 
CALL LAMNAR 
IF .( ERROR ) GO TO 10 
IF (SEPRM GO TO 20 
IF (.NOT. TRANS) GO TO 20 
CALL TURSLN 
IF (ERROR) GO TO 10 
520 CALL PROFIL 

ALPHl * 3.1A159/2.0 - ALPl 

CALL AFM1X( ALPHl, TE,SP,ME{NST),NP) 

CALL N0Z2LC(NSP, ALPHl, NP) 

GC TO 10 
END 
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SIKFTX: INPU 


DECK 


SUBRCUTIKE INPUT 

CCPPCN/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEN,KSMTH, 
lKSPLNrKLE,KATCH,CTHET,DLAM,TLAM,DTUR8,TTURB,KPRE,K6RAD»KS0E,KLAM, 
2K«AIN»KPROF,XtlOO),Yt 100),PRESI 100 } ,UE ( 100 ) ,MEf 100 ) ,POPTZ ( 100 ) , 
3VCVCR<10C ),TWAU 100) 

COPMCN/Cq, /ERROR, TRANS, SEPRN 
LOGICAL ERROR,TRANS,SEPRN 
REAL PE 

ERROR* .FALSE. 

TRANS=; .FALSE. 

SEPRN* .FALSE. 

C7HET = C.O 
KATCH = C 
KSPLN = 1 
KLE = 0 
WRITEC6,1000) 

READ t5,1020) R,PTZ,TTZ,UPMACH 

WRITE(6,1050) 

MRITEt6,lC60) R , PTZ , TTZ , UPMACH 

READI5, 1C10)NVP,N7URB,KEM 

WRI7E86,1070)NS7,NVP,NTURB,KEM 

READ €5,1020) OL AW,TLAM,DTURB,TTURB 

WRI7EI6, 1080) DLAN,7LAM,07URB,7TURB 

READ 15,1010) KPRE,KGRAO,KSDE,KLAM,KMAIN,KPROF 

WR17e(6,109Q) KPRE,KGRAO,KSDE,KLAM,KMAIN,KPROF 

IF(NS7.G7. 100. OR. NTURB. G7.NST.0R.KEM.lt. 0. OR. KEM, 67. 1. OR. KSPLN.L7. 
lO.CR.KSPLN.GT.l.OR.KLE.LT.O.OR.KLE.GT. l.OR.KATCH.LT.O.OR.KATCH.GT. 
211 GC TO 70 
00 1 I »1,NST 
1 TWALU) * TTZ 

WRITE(6,1130) (X(I ),Y(I ),POPTZm,TWAL(I),I = l,NST) 

RETURN 

70 ERROR = .TRUE. 

MRITE(6,1I70) 

RETURN 

lOOC ECRKATI IFl////) 
iOlO FORMAT! 16 15) 

10S2C FORMAT! 8F 10. 5) 

1030 F0RMAT!4F10.5) 

1030 FORMAT! IFO, 6X,22HB0UN0ARY LAYER SECTION) 

1060 FORMAT!// 9X, IHR, IIX , 3HPTZ, 8X, 3HTTZ ,8X»6HUPMACH/ 3X 

1,F9.2,3X,F10.2,2X,F9.2,4X,F8.4) 

1070 FORMAT! /6X,3HNST,8X, 3HNVP,9X, 5HNTURB,7X,3HKEM/5X, 13, 8X, 13, lOX, 

I 13,9X,I2) 

1080 FORMAT! /6X,4HDLAM, 7X,4HTLAM,8X, 5HDTURB,7X,5HTTURB/4X,F10. 6,1 X,F10. 
16,2X,F10.6,2X^F10.6) 

1090 FORMAT! /6X,4HKPRE,7X, 5HKGRAD,7X,4HKSDE,8X,4HKLAM,7X,5HKMAIN,7X,5HK 
IPRCF/7X, I2,9X,I2, lOX, 12, lOX, 12, 9X, 1 2, lOX , 1 2 ) 

1130 FORMAT! /9X , IHX, 1 IX, IHY, lOX, 5HPQPTZ, 10X,4HTWAL/ ! 3X,F10.5 ,2X,F10.5 ,4 
lX,FlC.6,AX,F10.4n 

1170 FORMAT! ////10X,48HERROR IN INPUT DATA. RECHECK INPUT INSTRUCTIONS 

II 
END 
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SIBFTG PREGfl DECK 


SUBBCUTINE PRECAL 

CG(^MCN/Cl/GAM,R,PTZ,TTZ»UPHACH,NSTfNVP,NTUR8,KPVM»KEM,KSMTH, 
IKSPLN,KLE,KATCH, CTHET,DLAM,TLAM,DTURB,TTURB,KPREfKGRADtKSDE,KLAM, 
2KHAIN,KPR0F,X<100),YI 100),PRESI 100) jUEClOO) , ME (100) ,POPTZUOO) , 
3VCVCRU0C ) ,TWAL< 100) 

CCMMCN/C2/PSZ»TS2»UZ,ASZ,ATZ,RHSZ»RHTZtMUSZ»MUTZ,NtSZ,NUTZ»CPt 

1PR,TC,ARCL 

CQMMCN/C3/X0M(100),YOM(100)f S( 100),SOLU0O> ,AE( ICO) »TSE( 100) , 
ITAWLt 100) »TAWT« 100 ) ,TBAR i 100 ) ,RWi 100) , SW ( ICO) ,SUTHU100) , 
2RHSWU00),RHSE(100),HEADW( 100) fHEAOEClOO) ,NUW(100) t«UBAR( 100) , 
3AAU00) fEB< 100),FF( 100) , DUDS C 100), DMDSt 100) ,OMOH1CO) 

COMMON/ C9/ ERROR* TRANS tSEPRN 

DIMENSION SOERt 100) , CMU ( 20 ) , CPR I 20 ) ,CTC ( 20 ) 

REAL MUSZ*MUTZ,NUSZ,NUTZ,MUSLE,MUSLM,ME,NUW,MUBAR 
LOGICAL EPROR^TRANSjSEPRN 

READ DATA FCR MU* PR* AND TC CURVE FITS 

DATACCMUn ) ,I»1, 5)/-. 01945170, 1, 3019531, -.3A511323* 

1. 068277826, -.00566593/ 

DATAICPRU ) , I -1,5)/. 8557, -.2341 36*. 1078624, 

1-. 0236214, .00202863/ 

DATAICTCI I ) , I*il, 51/-. 0383932 3, 1.2697427, -.3091 1252* 

1.08743781, -.009674725/ 

INITIALIZE STATIC AND TOTAL PARAMETERS 

TSLE^ 51E.688 
TSLM= 288.160 
MUSLE- 3.711402E-7 
MUSLM- I.777029E-5 
TCSLE=? 3.202206E-3 
TCSLM:? 2.561796E-2 

TSZ=fc TTZ/I l.+(GAM-l. )/2.*UPMACH**2) 

PSZ* PTZ* i TSZ/TTZ)**! GAM/ ( GAM-1. ) ) 

RHSZ= PS2/R/TSZ 
8HTZ* PT2/R/TTZ 
ASZ5. SQRT(GAM#R*TSZ) 

ATZ= SQRTIGAM*R*TTZ) 

UZ= UPMACH*ASZ 
CP= R4GAM/IGAM-1.) 

IF IKEM.EQ.l) GO TO 10 
TCCN= 198.60 
TRl« TSZ/TSLE 
TR2» TTZ/TSLE 
GO TC 20 
10 TCCN^ lie. 33 
TR1=? TSZ/TSLM 
TR2= TTZ/TSLM 

20 CALL CURVFT(CPR*PR,TR1, 0,4,0) 

CALL CURVFT(CTC,TC,TR1, 0,4,0) 

CALL CURVFT(CMU,MUSZ,TR1, 0,4,0) 

CALL CURVFT(CMU,MUTZ,TR2, 0*4,0) 

IF (KEM.EQ.l) GO TO 30 
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TC= TC*TCSLE 
ML’SZ* MUSZ + MUSLE 
MUTZ^ MUTZ*MUSLE 
GG TC 40 
30 TC= TC*TCSLM 

MUSZ* MUSZ*MUSLM 
MUTZ= MUTZ*MUSLM 
.40 NUSZ= MU'Z/RHSZ 
NUTZa MUTZ/BHTZ 

C/^tCUL/STE GECyETRY RATIOS AND ARC LENGTHS 

XCKM)= >m/X(NST) 

Y0N(1)= Ym/X(NST) 

S11)*; 0. 

00 50 1=2, NST 
XCf'd)* Xm/XINST) 

YO>'< 1)= Yl I)/X(NST) 

50 S(I)=» S(I-U+SQRT( (Xm-X{ I-1J)**2+(Y{ n-YU-1) )**2) 

ARCL* S(NSTJ 
DC 60 1=1, NST 
*>0 SOU n= SI D/ARCL 

CALCULATE PRES ,UE,ME,POPTZ , AND VOVCR AT EACH STATION 

PRESSURE OVER TOTAL PRESSURE GIVEN AS INPUT 
13C DC 140 1*1, NST 

IF(PCPTZm.LT.O..OR.POPTZ(I).GT.l.»GO TO 290 
PRESm= POPTZm*PTZ 

UE( n=! SCRT(2.*GAW/{GAH-1.)*PTZ/RHTZ*I1.-<PRES(I)/PTZ)**( (GAM-l.I/ 
IGAPm 

TSEIII* ITZ>UEm**2/I2.*CP) 

AE(n*t SCRT(GAM*Rj»TSEU ) ) 

KE( UEI n/AE< I) 

I4C VOVCR 1 1 )= SQRTUGAM+1.»/(GAM-1.»*(1.-IPRES(I)/PTZ»**( IGAM-1. )/GAM) 

*n 

PRINT INITIAL CALCULATED PARAMETERS 
ITO WR1TE(6,1000) 

WRITEf6,1010> PSZ,^TSZ,.UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP, 
1PR,TC,ARCL 

WRITE (6,1020) II ,PRESm,UEIl),ME(I),POPTZ(n,VOVCRm,l = l,NST) 


PRINT GEOMETRY PARAMETERS 
200 IF IKPRE.NE.l) GO TO 210 

WRITE (6, 1030) I I,X( I) ,Y( I),S(I),XOM< I ) ,Y0M( I) ,SOL(I ) ,I=1,NST) 

CALCULATE OTHER NECESSARY PARAMETERS AT EACH STATION 

210 DO 220 1=1, NST 

TEMl^ 1.*.5*{QAM-1.)*ME( I)**2 
RHSWII)= PRESm/R/TWAL( I) 

RHSEU)= PRES(I)/R/TSE( I ) 

HEADWm* .5*RHSWII)+UEm**2 
HEADEin* .5*RHSEm*UE( I)**2 
SWI I )^ TWALI I )/TTZ-l. 
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SUTHUI SQRT(TWAUn/TTZ)*(TTZ+TCONI/(TWAL( n + TCONI 
NUW( I)« SUTHL (I)*NUTZ*{ l. + SWH ) )+*2*TEMl*» IGAM/ (GAR-1. » ) 

Rwm= UE( i)*sn >/Nuwm 

TAWLI n= TSEU)»( l.+PR**(l./2.)*(TEMl-l. ) J 
TAb'TC II* TSE( I)*(l.+PR**(l./3. )*(TEMl-l.l ) 

T8ARIII* .5*{TWALm+TSEn)l+,22*PR*#(l./3.}*ITTZ-TSE(IH 
RUBARt! )« RUTZ*SUTHL( I)*TBAR( II/TTZ 

iB( II* ME( I 1*ATZ/NUTZ*(TSE(I l/TTZI**((GAM+l.l/(2.*GAM-2.n 
AA(i )= een i*Tsen i/tbar( i i*(mubaru i/mutzi**.268 

FR(I }*1.*.1599*MEMI**2 + .60*SW( I I + .2101+SW ( I I*ME (I I **2 + .0114*ME 1 1 1 
I**A*.018C*SW( II*ME( Il**4 + .l825*SHni**2+.0735*SW(II**2*ME(I l**2 
2+.CQ73*Sh( I I**2*ME( 1 1**4 

2520 CONTINUE 
C 

C GORPUTE VELCCITY AND MACH NUMBER GRADIENTS ALONG THE SURFACE 
C 

DUCSINSTI * 0.0 
OMCSINSTI * 0.0 
I * 1 

2521 IF (I.EQ.NST) GO TO 240 
IP(MEUI .NE. ME(I + lll GO TO 222 
DUDSm * 0.0 

oMcsni * 0.0 
I = 1+1 
GO TO 221 

322 II* I 

323 I * I+l 

IFd.EQ.NSTl GO TO 224 
IF CRE( II.NE.MEI I+ll I GO TO 223 
224 NSP * I+l -II 
C SPLINE CURVE TECHNIQUE 

230 CALL SPLINE(StIlUUE(IllfNSPtOUOS<Il}»SOER(IUI 
CALL SPLINE(S(I1 1,ME( I1I,NSP,DM0S( 111, SDERdll) 

IF < I.NE.NSTI GO TO 221 
240 DC 250 1*1, NST 
250 DMDLdl* ARCL*DMDSd I 
C 

C PRINT OTHER CALCULATED PARAMETERS 
C 

IF(KPRE.NE.ll GO TO 260 

WRITE (6, 10501 d ,AE( 1 1 , TSE ( 1 1 , TWAL d I , TAWL d I ,TAWTd 1 ,TBARd 1 , 

11*1, NSTl 

WRITE (6, 1060 1 d ,RWd>,SWI 1 1 ,SUTHLd I ,RHSWd I ,RHSEd 1 ,HEADWd I , 
IHBAOEU l,NUW( II#MUBAR( I ), 1*1, NSTl 
260 IFiKGRAD.NE.ll GO TO 270 
WRITEI6, 10701 

WRITEf6,10801 (I ,DUDSd l,DMDSd IfDMDLd ),1*1,NSTI 
C 

C CHECK FOR IFPROPER INPUT 
G 

270 DC 280 1^2, NST 

IF 1UE( I l.NE.O. J GO TO 280 
ERROR’S .TRUE. 

MRITE(6, 10901 
RETURN 

280 CCNTINUE 
RETURN 

290 ERROR* .TRUE. 

MRITEt6,lI00l 

RETURN 
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c 

C fCRf«AT STATEMENTS 
C 

1(MJC FGRMAT< ll-l///4X,24HPRELIMINARY CALCULATIONS///) 

lOIC FORMAT! SX^IOHPSZ = F12.5/5X, lOHTSZ * F10.4/5X, lOHUZ 

I Fll.5//5X,10HASZ = F11.4/5X, lOHATZ = Fll.4//5X,10HRHSZ 

2?= G15.7/5X, lOHRHTZ = G15.7//5X,10HMUSZ » G15.7/5X, lOHMUTZ 

3 = G15.7//5X,10HNUSZ = G15.7/5X, lOHNUTZ = G15.7//5X, lOHCP 

A F11.A/5X,10HPR = F9.5/5X, lOHTC = G15.7/5X,10HARCL 

5 * P8.A///) 

I®20 FORMAT! /IX, 7HSTATION,7X, AHPRES,13X,2HUE, 12Xf2HME,l IX, 5HP0PTZ,9X,5H 
IVQVCR/! 2X, I3,5X,F12.5,3X,F12.5,AX,F10.6,AX,F10.6,AX,F10.6)) 

1030 FORMAT! /// IX, 7HS TAT ION, 7X, IHX, 12X, IHY, 12X, IHS,12X,3HX0M,9X,3HY0M, 
19X,3HS0L/!2X, 13, 3X,F12.5,1X,F12.5,1X,F12.5,4X,F9.5,3X,F9.5,3X,F9.5 
2)1 

1050 FORMAT! ///IX, 7HSTAT ION, 5X, 2HAE , lOX , 3HTSE , 9X ,4HTWAL , 8X ,4HTAWL , 8X ,4H 
1TAWT,8X,4HT8AR/! 2X, I 3 , 4X , F9. 3, 5 ! 4X ,F8 . 3 ) ) ) 

1060 FCRMAT!///1X,7HSTATI0N, 11X,2HRW,6X,2HSW,4X,5HSUTHL,7X,4HRHSW,12X,4 
IHRHSE,8X,5HHEADW,4X,5HHEADE,9X,3HNUW,12X,5HMUBAR/!2X,I3,3X,F16.1,2 
2X»F4.1, 1X,F7.3,2X»G14.6,2X,G14.6,1X,F8.3,1X,F8.3,2X,G14.6,2X,G14.6 
3) ) 

1070 FORMAT! 1K1///21X,17HSURFACE GRADIENTS///) 

1080 FORMAT! IX , 7HSTAT ION , 13X , 4HDUDS, 15X ,4HDMDS , 15X ,4H0MDL/ ! 2X, 13 ,4X,F18 
1.6,IX^F18.6,1X,F18.6)) 

1090 FORMAT! /////lOX, 83HTHERE IS A STAGNATION POINT AT A STATION OTHER 
ITHAN STATION 1. THIS IS NOT ALLOWED) 

HOC FORMAT! /////lOX, lllHAN INPUT PRESSURE, VELOCITY, OR MACH NUMBER IS 
lEITHER LESS THAN ZERO OR GREATER THAN ITS MAXIMUM ALLOWABLE VALUE) 
ENC 


$IBFTC LAMNA DECK 

SUBRCUTINE LAMNAR 

COMMON/ C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH, 
1K3PLN,KLE,KATCH,CTHET,DLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM, 
2KMAIN,KPRQF,X1100),Y! 100), PRES! 100 ) , UE ! 100 ) ,ME ! 100 ) ,P0PTZ!100) , 
3VOVCRnOC),TWAL! 100) 

C0MMCN/C2/PSZ,TSZ,UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP, 

1PR,TC,ARCL 

CfMMCN/C3/XOM!100),YOM! 100 ) , S! 100 ) , SOL ! 100 ) , AE ! ICO) ,TSE!100) , 
ITAWL! 100),TAWT( 100) ,TBAR!100),RW!100),SW11C0) ,SUTHL(100) , 
2RHSW!100),RHSE! 100),HEADW! 100 ) ,HEADE ! 100 ) ,NUW( 100) , MUBAR ! 100 ) , 
3AA! ICO) ,EB! 100), FF ! 100) , DUOS ( 100 ), DMDS ! 100 ) ,DMDL ! 100 ) 
COMMON/C4/THET!100),OELSR! 100) , DELTA! 100) , FORM! ICO) , 

IFQRMHIOC ) ,FORMTR( 100 ) , RTH! ICO) ,RTHI ! 100) ,CF!100) , 

ITAUW! 100),NUSS! 100 ) , OTDY ! 100 ) ,HTRAN ( 100 ) ,CRN ! 100 ) 

CCMMGN/C5/SHAPLC 100 ) , SHAPK ! ICO ) ,8,NS 
CGMM0N/C6/FTRAN, FORMS 
COMMCN/CT/INST, ITRAN,ISEP 
COMMON/CS /ERROR, TRANS, SEPRN 

DIMENSION CORLNI 100 ), CORML ! 100 ), SHEAR! 100 ) ,DTH( 100) 

DIMENSION CCN!20),CRCR! 20),CDIF!20),CSHR(20),CCRN!20) ,CDTH!20) 
DIMENSION STAB!505),CTABl!505),CTAB2!505) 

REAL MUSZ,NUSZ,MUTZ,NUTZ,ME,NUW,MUBAR,NUSS,NURW,K8AR,INT1,INT2 
LOGICAL ERROR,TRANS,SEPRN 
EXTERNAL FUNCT, INTI, INT2 
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c 

C HE/II]) DATA FCR CORLNU), RCRIT DIFF, SHEAR, CRN, AND DTH CURVE FITS 
C 

DATAICCNI I J,I=?1,6»/-.08178,.C6670,-.03143, 

I. 00873, .C 1657,-. 01052/ 

DAT A ( CRCP ( n, 1=1,6 )/5, ^7073, A3. 6053, 227. 198, 
1-2067.04,-27172.7,13691.2/ 

BATAICDIFi I), 1=1, 6)/903. 785, 26365. 0,3. 85695E+5, 

II. 110A4E46,-A.53853E+7,-7.70276E+7/ 

DATAICSHRC I) , I=l , 16 )/. 224488, -1.91539, -9. 894, -68. 13488, 

1-. 00 151 2, -1.4768, -10. 52925, -152. 2781, -.002406, -.015629, 
1-1.45743,-126.23395, .000752, .005385,. 917838, -39. 40644/ 

DATAICCRM 1), 1=1, 161/2.02056, -19. 721 1,-24. 0495, -1400. 002, 

I-. 050979, -10. 880 12, 62. 4419, -508 1.76, -.014343, 2. 279845, 
1129.7008,-6257.848, .0270567,-1.677051,57.4397,-2552.266/ 

OATA(COTH I), 1=1, 161/8.02829,-4. 30978, 88. 8244, 36. 4336, 
12.71101,-7.42259,242.293,-16.293, -.16394,-7.61942,286.9795, 

164. 11 186, -.16758, -3. 70289, 130. 8107, 11 1.3276/ 

C 

C INITIALIZE PARAMETERS 
C 

INST = 0 
ITRAN = C 
ISEP =? 0 
CFU)^ 0. 

TAUWm* 0. 

NUSS<1)= 0. 

DTCY(1)= 0. 

MTRAMDs 0. 

CRN(1)= C. 

RTRANs; 0. 

RLE » 0 
C 

C CHECK CONSISTENCY OF INITIAL VALUES 
C 

IE (DLAM.GE.O..ANO.TLAM.GE.O..ANO.DTURB.GE.O..AND.TTURB.GE.O.) 

IGC TO 10 
ERROR = .TRUE. 

WBITEt6, 1000) 

RETURN 

10 IF (NTURE.NE.l) GO TO 30 
ITRAN = 1 

If (DTURE.GT.O.. AND.TTURB.GT.O. » GO TO 20 
ERROR = .TRUE. 

MRITEL6,1010) 

RETURN 

20 IF (Uem.GT.O.) GO TO 240 
ERROR = .TRUE. 

WRIT|{6, 1020) 

RETURN 

C 

C BEGIN CALCULATION IN LAMINAR REGION - CHECK FOR INITIAL VALUES 
C CALCULATE INITIAL CORRELATION NUMBER 
C 

30 IE (CLAM.EC.O..AND.TLAM.EQ.O.) GO TO 70 
IF (Ufl D.GT.O.) GO TO 40 
ERROR = .TRUE. 

WRITEI6,1030) 

RETURN 

40 IF ITLAM.EC.O.) GO TO 50 
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C INITIAL MOMENTUM THICKNESS WAS GIVEN 
T6M1= l.*.5*(GAM-l. J*MEm**2 

CGRMUl )* -ATZ*TLAM**2/NUTZ/SUTHL( I ) /ARCL/TEMl** (I3.-GAM) / 
l(2.*GAM-2.n 

CORLMl) » C0RML(1)*0MDL(1J 
G13 TC 90 

C INITIAL DISPLACEMENT THICKNESS WAS GIVEN 
SC IF (ABSICMDLIin.GE. .00011 60 TO 60 
C0BLN«1)= 0. 

TEMl* 1.4.5*(GAM-I.)*ME( l)»*2 

FORMID* 2.38All+U. + <2.79-1.78*PR**.5)*m. + SWIl))*TEMl-l.) » + (4.6 
r5*PR**a./3.)-3.65*PR».5)*PR**.5*lTEMl-l.) 

THETIU* OLAH/FORMm 

CORMLd )» -ATZ*THET( l)**2/NUTZ/SUTHLm/ARCL/TEMl**n3.-GAM)/(2.*G 
I AM-2.)) 

GO TC 90 

60 IF (CMDUD.GT.O.) CALL ROOT (-1. ,0. «DLAM,FUNCT, . 5E-5 1 CORLNI I ) ,SL ) 
IF (EMDL(l).LT.O.l CALL ROOT! 0. , . 2.DLAM ,FUNCT, . 5E-5 tCORLN II ) , SL) 
CQRMLIl ) ^ CORLNI l)/OMOU 1) 

GO TC 90 

NC INITIAL LAMINAR VALUES GIVEN 
CALCULATE INITIAL CORRELATION NUMBER 

SFARP LEADING EDGE 

7C IFIKLE.NE.l. AND. A6S<DMDL(1)).GE.. 0001) GO TO 80 
CGRLN(1)= 0. 

C'QRML(l)* 0. 

GO TC 90 

STAGNATION POINT 

8C CALL CURVFTICCN,CORLN(l),SWm, 0,5,0) 

COBMU l )* CORLNi 1)/0MDLU) 

IF ICCRMLI D.LT.O. ) GO TO 90 
ERRCR=? .TRUE. 

WRITE(6,10A0) 

RETURN 

SOLVE LAMINAR DIFFERENTIAL EQUATION 
CALCULATE CORRELATION NUMBERS ALONG THE SURFACE 

90 TBMl* l.T.5*(GAM-l.)*ME( 1)**2 
TEM2= (3.*GAM-1. )/(2.*GAM-2.) 
iBL» 0.0C2*ARCL 
SS* -CEL 
NTAE=1 

CTAEltl)= CORLNI 11 
CTAB2I.1 )= CORMLI U 
STAB(1)= 0. 
too SS» SS+DEL 

SSDEL = SS+DEL 

CALL LGRNQE(S,SW,NST,SS,ANS1) 

CALL LGRNCE(S,ME,NST,SS,ANS2) 

CALL LGRNGEIS, ME, NST, SSDEL, ANS3) 

CALL LGRNGE(S,DMOL, NST, SSDEL, ANSA) 

Al= 0.4363 1-0. 00367* ANSl+0.00A81*ANSl**2+0. 0065 1*ANS1**3 
A2* 5.43220+2.25400*ANSl-0.06672*ANSl**2-0.20637*ANSl**3 
A3» 4.51903-10.49775*ANS1-12.71732*ANS1**2“2.95270*ANS1**3 
A4* 19.01831+62.76597*ANS1+115.00986*ANS1**2+62.53113*ANS1**3 
A=^ A1-A3*CTAB1(NTAB)**2-2.*A4*CTAB1INTA8)**3 
B» A2 + 2.*A3*CTABHNTAB)+3.*A4*CTAB1(NTAB)**2 
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ROR SW — C mO 

IF { CTAEUNTAB) .CE.-.l) GO TO 101 
A^.3953 
8^4.739 
101 K1 = 0 

SOLI H SS/ARCL 
S0L2 =5 SSDEL/ARCL 
TiK3 =i SIMPSliSOLI,SOL2, INT1,K1) 

IF {TEM3.EG.O..OR.K1.EQ.O) GO TO 110 
ERRCR=! .TRUE. 

WRITE(6,1050) 

RETURN 

lie IF INTAB. GT.l) TE>1A= ANS2**( -B )*TEM1**TEM2 
TBMl® 1.4.5*(6AM-1. )+ANS3**2 
TBN5= ANS3**(-8)*TEM1**TEM2 
TER6* -A*TEN5*TEM3 
IF (NTAB.EC.l) TEW7=C. 

IF (NTAB.GT.l) TEW7= TEM5/TEM4*CTAB2(NTAB) 

NTABs NTAB+1 
CTAE2(NTAB)= TEM6+TEM7 
CTAB1{NTAB)= CTAB2 I NTAB ) *ANS4 
STABINTAE)= SSDEL 

WHEN SW IS NOT EQUAL TO 0.0 tCURVE FIT RANGE ON CORLN 
IS FRCM -.32 TO .16 

IE (CTABIINTAB).GT. .16) GO TO 120 
IR (SS.LT.ARCL) GO TO 100 
1220 IF (KSOE.NE.l) GO TO 130 
WRITE66,1C60) 

WRITE (6, 1070) (STAB! II, CTAB III), 1=1, NTAB) 

CALCULATE LAMINAR BOUNDARY LAYER PARAMETERS AT EACH STATION 

130 IF (KLAM.NE.l) GO TO 140 
WRITEI6, 1080) 

E40 1= 0 
150 1= 141 

IF (I.EG.NTURB) ITRAN»-1 
IF (Sfl ) .LE.STAB(NTAB) ) GO TO 160 
ERROR= .TRUE. 

WRlTE(i5,1090) 

RETURN 

r«0 CALL LGRNGEISTA8,CTABl,NTAB,Sm,C0RLN(m 
CALL LGRNGEI STAB, CTAB2, NTAB, SI I),CORMH I) ) 

•CBIAIN SHEAR, CRN, AND DTH FROM CURVE FITS VS CORLN AND SW 
CALL CURVFTICSHR, SHEAR! I ),CORLN( I) ,SW1 I ) ,3,3) 

CALL CURVFT(CCRN,CRN( I) ,CORLN( I ),SW( I ),3,3) 

CALL CURVFT(CDTH,DTH(I),CORLN(I),SW(I),3,3) 

FOR SW = C.O 

IRICORLNII ).GE. -.1) GO TO 161 
SHEAR(I) =1 -1.2222*C0RLN( I ) + .26 
CRNtl) = -58.824 #CORLN ( I )-.6765 
DTHII) = -22.222*CORLNI I )+7.1112 
C CALCULATE OTHER LAMINAR BOUNDARY LAYER PARAMETERS 
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lai TEMl» l.-*.5*<GAM-l.y#MEU)**2 

THET( I)= SQRT(-CORML( I )*NUTZ*SUTHL ( I )*ARCL/ATZ*TEy«l**C (3. -GAM)/ 
I(2.*GAM-2. ))) 

FORM! I)= C-l.ll38=*C0RLN! I) +2 .38411 )*( 1 . + ( 2.79-1 . 78*PR**. 5 ) ♦< (1.+ 
ISW! I ) )*TEM1-1.) )+(4.65*PR**! 1./3. )-3.65*PR**.5)*PR**.5*(TEHl-l. ) 
DELSR(U» THETCI )4FQRM( I) 

RTH(I)= ten )*THET! n/NJUW( I ) 

FCRMId )s (FORM! I )-SQRT ( PR )♦ ( TEMl-1. ) ) / ( ( I .+SW (I ) )*TEMi) 

FORMTR! F0RMI(I)*(1.+SW( I)) 

DELTA! I )= THET! I )*1DTH( I )+(TEMl-l. )* IFDRMTR ! I)+l . ) ) 

SHAPLCI )» DELTA! I)**2/NUW! I)*DUOS! I ) 

IF II.EC.l) GO TO 180 

CFRW* 2. ♦SPEAR! I )*SQRT! -SOL ! I) /ME! I ) /CORML ! 1) ) 

CR!I»* CFRW/SQRT!RW! n > 

TAUW!I)= CF( I)*HEADW! I) 

NURW^ CFPW*PR**.3/CRN( n 
NUSSm* NURW*SORT!RH(I) ) 

DTCYCn* NUSSin*!TAWL( n-TWAL! in/su) 

HTRAMII* TC*DTDY!I) 

IF (TAUWin.GT.O.) GO TO 180 
IF IKATCF.NE.O) GO TO 170 
ISEP? I 
SBPRN= .TRUE. 

RETURN 

170 ITRAN* -2 
GO TO 2 70 

18C IF ! I.EQ.l.AND.UE! D.EQ.O. ) GO TO 190 

SHAPKd )* NUTZ*RTH( I )**2*SUTHL! I )**2*! l.+SW d ) ) ♦♦4/ATZ/ME ! I) **2/ 
IFF !I)/ARCL*DMDL! I)*TEM1**! I./IGAM-I. ) ) 

GO TO 20C 

190 SHAPKd)* 0.07 

200 RTHim* RTH(I)*SUTHLm*{l.+SWin)**2/FFd)/SQRT(TEMl) 

CALCULATE RCRIT TO CHECK FOR INSTABILITY AND TRANSITION 

CALL CURVFT !CRCR,RCR IT, SHAPK ! I), 0,5,0) 

IF ISHAPKtl) .GT. .07) RCRIT * 8.3163 
RCRIT* EXPIRCRIT) 

IF! INST.NE.O) GO TO 210 

CHECK FOR INSTABILITY 

IF!RTHH II.LT.RCRIT) GO TO 270 
RINS= RTFI! I ) 

INST* I 
GO TC 270 

CHECK FCR TRANSITION 

210 Kl= 0 
NS= I 

TEM= SI MP3 1! SOL I INST), SOL! I),INT2,K1) 

IF ITEM.EQ.O..OR.K1.EQ.O) GO TO 220 
ERROR* .TRUE. 

MRITEt6,1100) 

RETURN 
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?20 KBAR=? TE(»/ (SOLU J-SOLIINST) ) 

CALL CURVPHCDIF,OIFF,KBAR,0,5,0» 

IR fKBAR .GT. .031 DIFF =* 44COO.*KBAR+700.0 
R1RAN=! RINS+DIFF 
IRIRTHI { I i.LT.RTRAN) GO TO 270 
IF U.LT.NTURB) GO TO 270 
ITRANa -1 
GO TO 27C 
23C ITRAN=! I 

GOfsPUTE INITIAL VALUES FOR TURBULENT SOLUTION 
240 TRANS= .TRUE. 

IF ICTURE.EQ.O.. AND.TTURB.EQ.O. ) GO TO 260 
IF {CTURB.GT.O.. AND.TTURB.GT.O. ) GO TO 250 
ERROR = .TRUE. 

WRITEC6,1110) 

RETURN 

250 THETUTRAN)= TTURB 

FORN{ITRAN)= DTURB/TTURB 

TERl=e l.*.5*IGAM-l. )*ME( ITRAN)**2 

FCRNH ITRAN)* (FORM! ITRAN )-PR** ( 1. /3. )*(TEM1-1 . ) )/( ( l. + SWI I TRAN) ) 
1*TEPU 

260 IF (CTNET.GT.O.. AND. DTURB.EQ.O.. AND. TTURB. EQ.O. ) THET(ITRAN)= 
1CTHET*THET{ITRAN ) 

TMETTR= TMET( ITR AN ) * ( TSE (ITR AN ) /TTZ)** ( (GAM+l , ) / ( 2. *GAM-2 . ) ) 

FTRAN= (NEUTRAN )*ATZ*THETTR/NUTZ)**1.268 
IF (RTRAN.LE.O. ) RTRAN=1000. 

FORKS’? FCRKK ITR AN )-0. 59389-0. 06591* AL0G(RTRAN)+0.001272*(AL0G«RTR 
IAN))**2 

IF (CTURE.GT.O.. AND.TTURB.GT.O.) FORMS=FORMI (1 TRAN ) 

RETURN 

PRINT OUTPUT 

270 IF (KLAK.NE.l) GO TO 280 

IFIINST.EQ.O .OR. INST.EQ.I) WRITE(6,1120) I ,CORLN(I ) ,SHEAR( I ) , 
iOTHI I ),FCRMTR)I) ,SHAPL( I ),RTHI ( I ),SHAPK( I ) ,RCRIT 
IPdNST.NE.O .AND. INST.NE.I) WR ITE ( 6, 1130 ) I ,CORLN ( I ) , SHEAR (I ) , 
IDTH( I ),FCRKTR( I) ,SHAPL( I ) , RTF I( I ) ,KBAR,DI FF ,RTRAN 
IF CITRAN.EC.-2) MR ITEI 6 , 1 140 ) 

280 IF<ITRAN.EG.-1.0R.ITRAN.EQ.-2) GO TO 230 
IF n.EG.NST) RETURN 
GO TO 15C 


FORMAT STATEMENTS 

l&OO FORMAT! ///«//, 10X,60HA NEGATIVE INITIAL VALUE HAS BEEN GIVEN. THIS 
IIS NOT ALLOWED) 

rOIO FORMAT!/////, 10X,75HINITIAL VALUES WERE NOT GIVEN FOR THE TURBULEN 
IT BOUNDARY LAYER AT STATION 1) 

1020 FORMAT! /////, 10X,«0HINITIAL VALUES WERE GIVEN FOR THE TURBULENT BO 
lUNCARY LAYER AT A STAGNATION POINT) 

E030 FORMAT! /////, 10X,94HINITIAL VALUES OTHER THAN ZERO WERE GIVEN FOR 
ITHE LAMINAR BOUNDARY LAYER AT A STAGNATION POINT) 
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I0<f0 FORMATI /////, 10X,106HF0R THIS INPUT DATA STATION 1 IS ASSUMED TO B 
IE A STAGKATION POINT, SINCE NO INITIAL THICKNESSES ARE GIVEN./ 
210X,118HIN THIS CASE PRESSURE SHOULD DECREASE INITIALLY. EITHER G 
3IVE AN INITIAL VALUE FOR DISPLACEMENT OR MOMENTUM THICKNESS,/ 
A10X,60H0R BEGIN WITH A SHORT REGION OF FAVORABLE PRESSURE GRADIENT 
5.1 

1050 FORMAT! /////, IOX,37HERROR IN COMPUTING INTEGRAL FOR CORLN) 
ia60 F0RMAT(1H///7X, 50HLAMINAR DIFFERENTIAL EQUATION - SOLUTION FOR CO 
IRLN///5124H S CORLN )//) 

1070 FGRMATn5!F12.5,2X,F7.4,3Xn ) 

1080 FORMAT! 1HI///1X, 59HLAMINAR CALCULATION OF INSTABILITY AND TRANSITI 
ION LOCATIONS/// I X,7HSTAT ION, 2X,5HC0RLN,5X,5HSHEAR, 5 X,3H0TH,6X,6HF0 
2RMTR,4X,5HSHAPL, 9X,4HRTHI,6X,5HSHAPK,9X,5HRCRIT,6X,4HKBAR,10X,4HDI 
3FF,9X,5HPTRAN) 

1090 FORMAT! /////, 10X,65HLAMINAR SOLUTION HAS PROCEEDED BEYOND THE RANG 
IE WHERE IT IS VALID) 

riOO FORMAT!/////, 10X,36HERR0R IN COMPUTING INTEGRAL FOR KBAR) 
llIO FORMAT! /////, 10X.64HIF INITIAL TURBULENT VALUES ARE GIVEN, THEY BO 
ITH MUST EE NONZERO) 

1120 FORMAT! 14,1X,5F10.4,1X,F12.1,IX,F10.5,1X,F12.1) 

1130 FORMAT! 14 , IX , 5F10.4, IX, F 12. 1 ,24X,F 12. 5 , 1X,F12. 1 , 1X,F12 . 1 ) 

1140 FORMAT!/////, 10X,85HLAMINAR SEPARATION HAS OCCURRED. ASSUMED TO BE 
1 TRANSITION TO TURBULENT BOUNDARY LAYER) 

ENC 


SIBFTC TURBL DECK 

SUERCUTIKE TURBLN 

CGMMCN/C1/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTUR8,KPVM,KEM,KSMTH, 
lK3PLN,KLE,KATCH,CTHETrDLAM,TLAM,DTURB,TTURB,KPRE,KGRAD,KSDE,KLAM, 
2KMAIN,KPR0F,X!100),Y! 100),PRES!100),UE!100) ,ME!100) ,P0PTZ(100) , 
3VOVCRI10C ) ,TWAL! 100) 

C0MMGN/C2/PSZ,TSZrUZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ,MUTZ,NUSZ,NUTZ,CP, 

1PR,TC,ARCL 

COMM0N/C3/X0M!100),YOM! 100 ) , S! 100 ) , SOL ! 100 ) ,AE ! 100) ,TSE!100) , 
ITAWU100),TAWT!100),TBAR! 100),RW(100)»SW!100) ,SUTHL(100) , 

2RHSW! 100 ) , RHSE! 100 ) ,HEADW! 100) ,HEADE ! 100 ) ,NUW! 100) ,MUBAR! 100) , 
3AAI100) ,EB 1100), FF! ICO), DUDS (100 ),DMDS! 100 ),OMDL! 100) 
COMMCN/C4/THET!100),DELSR! 100) , DELTA! 100), FORM! 100) , 

IFQRMUilOC ) ,FORMTR! 100 ) , RTH! 100 ) , RTHI ! 100 ) ,CF ! 100 ) , 

ITAUW ! 100 ) , NUSS! 100 ) , DTDY ! 100 ) , HTRAN 1 100 ) ,CRN ! 100 ) 

GCMMON/Ce/FTRAN, FORMS 
CCMM0N/C7/INST, ITRAN, ISEP 

C0MMCN/CE/XTAB!50.5),YTABII505),YTAB2!505) ,NTAB 
COMMCN/C9/ERROR,TRANS,.SEPRN 
REAL MUSZ,NUSZ,MUTZ,NUTZ,ME,NUW,MU8AR,NUSS 
LOGICAL ERROR,TRANS,SEPRN 
C 

C SCLVe TURBULENT BOUNDARY LAYER DIFFERENTIAL EQUATIONS 
C USING RUNGA-KUTTA 
CALL RUNKUT 

IF !KSDE.NE.l) GO TO 10 
NRITE!6, 1000) 

WRITE!6, 1010) !XTAB! Ii,YTABlU ),YTAB2U ),I»1,NTAB) 
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CALCULATE TURBULENT BOUNDARY LAYER PARAMETERS AT EACH STATION 

10 DC 30 I=ITRAN,NST 

IF (Sn ).LE.XTAB(NTAB)) GO TO 20 
ISEP =? I-l 
SEPRN=! .TRUE. 

RETURN 

20 TEMl » l.-*.5*(GAM-l. )*ME( I )**2 

CALL LGRNGE(XTAB,YTABlfNTAB,S(I ) ,F> 

THETTR= NUTZ*F**.7886/ME( I I/ATZ 

THETI n* THETTR*(TTZ/TSE(I ) )♦♦( CGAM+1. )/(2.*GAM-2. )) 

RTHin* LE( n*THET( I)/NUW( 1) 

CALL LGRNGEIXTAB,YTAB2,NTAB,S< 1 1 ,FORMI ( I ) ) 

FeRMTR(n=» FORMi ( I )*(i.+swn n 

FORM! I) = FORMTRI n*TEMl + PR** ( 1 . /3. )♦ ( TEMl-1 , ) 

OBLSRm=‘ THET( n*FORM( I ) 

P01«ER=» 2.0/{F0RMI( n-l.O) 

IF (FORMIC n.LT. 1.02) POMER=lOO. 

OELTA(I)^ ( l.+POWER)*OELSR( I ) 

CF( I)= 0.2A6*EXP(-1.561*F0RMI( I ) )*(UE(I )*THET( I) /NUTZ/I TEMl**{l./( 
IttAM-1.) ) n*#(-.268)*TSE( I ) /TEAR { I ) ♦ ( MUBAR ( I) /MUTZ ) *« ( .268) 

TAUWm= CF(n#HEADE(I) 

IF (I.EQ.l) GO TO 30 

MTRANU )a CF( I ) / 2. /PR**( 2. /3. )*RHSE ( I ) *UE ( I ) *CP* ( TA WT( I) -TWAL (I ) ) 
DTCYU)= HTRAN( I ) /TC 

NUSS( I)» S( I)*DTDY( I)/(TAWT( I)-TWAL( m 
CRN(I)= CPU )»RW( I)/NUSS( I ) 

30 CONTINUE 
RETURN 

1000 FGRMAT(lH///5Xf 62HTURBULENT DIFFERENTIAL EQUATIONS - SOLUTION FOR 
1 F AND F0RMI///4(31H S F FORMI )//) 

1010 FGRMAT(CA(F10.5t2X,F8.1,2X,F7.4,2Xm 
END 
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$1BJ=TG PRCFl DECK 


SUBROUTI^E PROFIL 

COr«HON/Cl/GAM,R,PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSMTH, 
lKSPLW»KLE,KATCH,CTHETtDLAM,TLAM,DTURBtTTURB,KPRE»KGRAD,KSOE»KLAM, 
2KFIAIN,KPPGF,XC100)»Y(100),PRES( ICO ) .UE ( 100 ) tME { 100 5 ,P0PTZJ100) t 
3VQVCR(10CJ ,TWALI 100) 

CGPPON/C3/X0M(100) fYOM( 100) , SUOO ) , SOL 1 100 ) , AE ( 100) ,TSE<100 ) , 
lTAMLClOO),TAWT(100),TBAR(100),RW(100),SwaCO)tSUTHL(100), 
2RMSW(100)rRHSE( 100)»HEADW( 100) t HEAOEI 100) »NUWU00) ,MUBAR(100) , 
3AAUC0) ,BB( 100)iFF( 100) ,DUDS1100),0MDSa00),DMDU100) 
CQMMON/CA/THETdOO) .DEL SRC 100), DELTA (100) tFORM( 100) , 
lFaRPniOO),FORMTRnOO),RTH(100),RTHI (lOO).CF(lOO), 
1TAUW(100),NUSS( 100),OTDY{100),HTRANl lOO).CRN(lOO) 

CGRKCN/C5/SHAPLI 100) ,SHAPK{ ICO) ,B,NS 
CQRKGN/CT/INST, ITRAN, I SEP 
RIAL ME.KGSS 

RRINl LCCATIONS OF INSTABILITY, TRANSITION, AND SEPARATION 

III (KPAIN.NE.l) GO TO 60 
WRITE(6, 1000) 

LFI INST.EQ.O) GO TO 10 
WRITE(6,1010) INST 
GO TO 20 
10 HRITEt6,1020) 

20 IR (ITRAA.LE.1) GO TO 30 
HRlTEr6,1030) ITRAN 
GO TO AO 
30 NRITE(6,10A0) 

40 IF(ISEP.EQ.O) GO TO 50 
WRITE(6,1050) ISEP 
60 TO 60 

50 NRITE(6,1060) 

PRINT LCCATIONS OF LAMINAR AND TURBULENT BOUNDARY LAYERS 
60 IBNC =1 ITRAN-1 

IE ( lEND.EQ.-l.OR.IEND.EQ.O) IENO=ISEP 
IF (lEND.EC.O) IENO=NST 
IF (KMAIN.NE.l) GO TO 70 
IF (ITRAN, EQ.l) NR ITE(6, 1070 ) 

IF (ITRAN.NE.l) WRITE(6, 1080 ) lEND 
IF (ITRAN.EC.O) WR ITE( 6 , 1090 ) 

IF (ITRAN.EQ.l) WRITE(6, 1100 ) ITRAN, lENO 
70 IF ( ITRAN. LE.l) GO TO 80 
lENC ^ ISEP 

IF(IEND.EC.O) IEND=NST 
IF (KKAIN.NE.l) GO TO 90 
MRlTEt6,llC0) ITRAN, IBND 

PRINT CALCULATED BOUNDARY LAYER PARAMETERS 

80 IF (KMAIN.NE.l) GO TO 90 
MRITE(6,1110) 

WRITE (6, 1120) (I,X(I),S(1),DELSR(I),THET(I ),DELTA(I),FORM(I), 
IFGRMI(I),I*=1,IEN0) 

NRITE(6,1130) 

WRITE (6, 1140) ( I,CF( I),TAUW( I),RTH(I »,DTOY(I),NUSS(I),HTRAN(I), 
1CRN( I), 1=1, lEND) 
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eCNPUTE BOUNDS ON VELOCITY PROFILES 

90 IF lKPR0F.NE.il RETURN 
WRITE (6, 1150 I 
IFUTRAN.NE.O) GO TO 100 
IL1=!‘ 2 
IL2= lENC 
ITl* 0 
IT2?= 0 
GG TC lie 

100 IL1= 2 

112= ITR/N-1 
IT1= ITRAN 
IT2^ lENC 

IP (ITl.EQ.ll IT1*2 

CALCULATE AND PRINT LAMINAR BOUNDARY LAYER VELOCITY PROFILES 

HO NVP1= MVP + 1 

IF (IL2.LT.ILI) GO TO lAO 
DO 130 I>ILltIL2 
WRITE(6«1160) I 
AAA* 2.+SHAPL(I)/6. 

BBB= *.5>*SHAPL( I ) 

CCC* -2.4.5*SHAPL( 1 1 
DDC* l.-SHAPL(II/6. 

DEL* CELTAII)/FLOAT(NVP) 

YP* -CEL 
DO 120 J*l,NVPl 
YP* YP+DEL 
ETA* YP/CELTA(I) 

YXPAX* YP/XINST) 

UUE* AAA*ETA+8BB*eTA**2+CCC*ETA**3+DDD*ETA**4 
U=< UUE*UE(H 

120 HRITE(6,1180) ETArYP,YXMAX,U,UUE 

130 CONTINUE 

CALCULATE AND PRINT TURBULENT BOUNDARY LAYER VELOCITY PROFILES 

lAC IFdTl.EC.O) RETURN 
DO 160 I^IT1,IT2 
POWER* DELTAdl/OELSRID-l. 

WRITE(6,1170) I, POWER 
DEL* OELTAd)/FLOAT(NVP) 

YP* -CEL 
DC 150 J*1,NVP1 
YP* YP+DEL 
ETA* YP/CELTACI) 

YXPAX* YF/X(NSTI 
UUE* ETA**d. /POWER) 

U* UUE+UEII) 

150 WRITE(6,1180) ETA^YP, YXMAX.U.UUE 

1160 CONTINUE 
RETURN 
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FCRMAT STATEMENTS 

IGOC F0RMAT( 1H///1X, 36HPRINCIPAL BOUNDARY LAYER INFORMATION///} 
lOIC FORMAT (/10X,31HINSTABILITY OCCURS AT STATION ,13) 

102C FORMAT I /lOX , 26H INSTABIL ITY DOES NOT OCCUR) 

1030 FORMAT ( / lOX , 30HTR ANS IT ION OCCURS AT STATION ,13) 

104C FORMAT ( / lOX, 25HTRANSITION DOES NOT OCCUR) 

1050 FORMAT ( / lOX, 30HSEPARAT ION OCCURS AT STATION ,13) 

IC%C FORMAT I /10X,25HSEPARATI0N DOES NOT OCCUR) 

1070 FORMAT ( / lOX , 37HL AMJNAR BOUNDARY LAYER DOES NOT OCCUR) 

1000 FORMAT I / lOX , A2HLAMINAR BOUNDARY LAYER - STATIONS 1 TO ,13) 

1090 FORMAT (/ lOX , 39H TURBULENT BOUNDARY LAYER DOES NOT OCCUR///) 

HOC FORMAT ( / lOX , 35HTURBULENT BOUNDARY LAYER - STATIONS, 2X, 

113, 6h TC ,13///) 

1110 FORMAT! /IX, 7HSTAT ION, 8X, IHX, 12X , IHS , 12X, 5HDELSR , lOX , AHTHET , 1 IX , 
15HDELTA,11X,4HFORM,10X,5HFORMI) 

1120 FORMAT! 2X, 13, 3X, 2F13. 6, F1A.6 , IX, FI A. 6, IX ,F1A. 6 , 1 X,2F14.4) 

1130 FORMAT! /// IX , 7HSTAT ION, 6X, 2HCF, 13X ,4HTAUW, 1 IX , 3HRTH , 14X,4HDTDY, 
113X,4HNUSS,10X,5HHTRAN, 12X,3HCRN) 

1140 FORMAT! 1 5 , F14.5, 2X , F14. 5 , IX, F12 . 1, 5X, F14.2 ,2X,F 14.2 ,1X, 
1F14.4,2X,F13.3) 

1150 FORMAT! 1F1///1X, 17HVEL0CITY PROFILES///) 

IF60 FORMAT ! / IX ,7HST AT ION , IX , 1 5, 2X, 7HPR0F I LE /3X ,7HY/DELTA ,9X , 
11HY,12X,6HY/XMAX, lOX.lHU, 12X,4HU/UE) 

1170 FORMAT ! / IX , 7HST AT ION , IX , 1 5, 2X, 7HPR0F ILE ,28X , 2HN=, IX ,F6.2/3X,7HY/D 
1ELTA,9X,1HY, 12X,6HY/XMAX,10X,1HU,12X,4HU/UE) 

1180 FORMAT! 1X,F8.4,2X,2G15.6,2X,F9.2,6X,F8.4) 

ENC 


IBFTG RUNKU DECK 

SUBRGUTINE RUNKUT 

RUNKUT SOLVES SIMULTANEOUS FIRST ORDER INITIAL VALUE 
ORDINARY DIFFERENTIAL EQUATIONS 

CGMMCN/C1/GAM,R, PTZ,TTZ,UPMACH,NST,NVP,NTURB,KPVM,KEM,KSHTH, 
IKSPLN,KLE,KATCH,CTHET,DLAM,TLAM,OTURB,TTURB,KPRE,KGRAD,KSOE,KLAM, 
2KMAIN,KPP0F,X!100),Y! 100 ) , PRES ! 100 ) , UE ! 100 ) , ME ! 100) , POPTZ ! 100 ) , 
3VOVCRUOC ) ,TWAL( 100) 

CCMM.CN/C3/XOM!100),YOM! 100 ) , SUOO) , SOL ! 100 ) , AE ! 100) ,TSE!100), 
ITAWL! 100) ,TAWT! 100 ) , TBAR ! 100 ) ,RW ! 100 ) , SW ! 100 ) ,SUTHL!100) , 
2RHSW!100),RHSE!100),HEADW!100),HEADE!100),NUW!100),HU8AR!100) , 
3AA! ICO) ,8B!100),FF! 100) , DUDS! 100), DMDS! 100) ,OMDL!1CO) 
CGMMCN/CA/FTRAN, FORMS 
C0MMGN/C7/INST,ITRAN,ISEP 

CGMMGN/C8/XTABI505) ,YTAB1! 505), YTAB2 1505), NTAB 
DIMENSION YY12),RY!2),YINC12),D0T!2),RUK12,4) 

DCU8LE PRECISION XX,RX,YY,RY,RUK,DEL,DOT, 
1TEM1,TEM2,TEM3,TE«4,TEM5,TEM6 
REAL ME,NUW,MUBAR 
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SET DEL SPACING AND STORE INITIAL VALUES 

DEL* 0.0C2*S(NST» 

10 YYm*FTRAN 
VYI2J* FCRKS 
XX* SCITPAN) 

NV*2 

NTA6 =? 1 
YTAB1(1)=^ YY(l) 

YTAe2tl)» YY(2) 

XTAem* XX 

SOLVE FOR YYfl) AND YY(2) AT NEXT XX INCREMENT 

SAVE PREVIOUS YYIII AND YY(2) 

20 DC 30 J»1,NV 
30 RYU I* YYU) 

RX* XX 

CALCULATE NEW YYIll AND YYI2J 
00 90 L*1,A 

C PUT DIFFERENTIAL EQUATIONS IN THE FORM OF 
C FIRST DERIVATIVE = REMAINDER OF EQUATION 
CALL LGRNGE(S^M£,NST,XX, ANSI) 

CALL LGRNGE{StSW,NST,XX»ANS2) 

CALL LGRNGE(StAA,NST,XX»ANS3) 

CALL LGRNGE<S,BBtNST,XX» ANSA) 

CALL LGR^'GE^S^DM0S,NST,XX,ANS5) 

CALL LGRNGE(S^TBAR,NST,XX,ANS6) 

TEMl* l.<( l.+ANS2)*YY(2) 

TEK2=5= .IZ3 + EXP(~1.561*YY<2) )*ANS3 

DCT( 1)= 1.268*(-YY< 1 ) /ANS1*ANS5*TEMI+TEM2 ) 

TEM3* YYI2)*(YYI 2) + 1. )**2* ( YYI 2 )-l . ) 

TEM4* l.*ANS2*IYY(2)*YY( 2)+4.*YY(2)-l. )/( ( YY(2)+l.)*(YY<2)+3.) ) 
TEM5* tYY{2)*YY< 2 J-1. )*YY< 2 ) /YY ( 1) ♦ ( . 123*EXP (-1 . 561 + YYI2 ) )*ANS3) 
TEMfi*;' (YYI 2 )*YY< 2)-l. )/YY( !)*♦( .7886)*<.011*lYY(2)*l. )*IYY(2)-1.) 
1#*2/Y¥(2)**2*TTZ/ANS6)*ANSA 
DOT (2)* -ANS5*.5/ANS1*TEM3*TEM4+TEM5-TEM6 
C APPLY THE RUNGA-KUTTA SCHEME 
BO AO J*1,NV 

AC RUKU»L) * OELAOOTU) 

GO TC (5Cf50f70,90)» L 
50 DO 60 J*1,NV 
'60 YYU)a RY( J)+RUKUfL)/2. 

XX* RXtOEL/2. 

SO TC 90 
70 DO 80 J = UNV 
8C YYU)* RYU)+RUKU»L) 

XX* RX+DEL 
90 CONTINUE 

C INCREMENT TFB DEPENDENT VARIABLES TO OBTAIN NEW YY(1) AND YY(2) 

OG ICO J*l,NV 

YINCU) * (RUKU,l)+2.*RUKU,2)+2.*RUKU,3)+RUKU»A))/6. 

100 YYU)* RYU)+YINCU) 

IF IYYI2).GT.2.8) RETURN 
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fSTGRB NEW CCMPUTED VALUES IN A TABLE 

NTAB ^ NTAB+1 
VTABltNTAB )= YY( 1) 

Y7AB21NT4BH YY(2i 

XTABINTABJa XX 

IF (XX.LT.SINST) ) GO TO 20 

RETURN 

END 


SIBFTC .FUNC DECK 

SUBROUTINE FUNCT ( XX , FX, DFX , INF ) 

GGFNGN/C1/GAM»R, PTZ,TTZ»UPMACH,NST,NVPtNTURB»KPVM,KEM,KSMTH, 
lKSPLN,KLE,KATCH»CTHET,0LAM,TLAM,0TUR8,TTURB,KPRE,KGRAD,KSDEtKLAM, 
2KMAIN»KPR0F,X(100IfY{ 100),PRES( 100»,UE(100),ME(100J .POPTZaOOI » 
3VGVCRI10C) »TWAL( 100> 

GGNWCN/C2/PSZ,TSZ»UZ,ASZ,ATZ,RHSZ,RHTZ,MUSZ»MUTZ,NUSZ,NUTZ»CP, 

1RR,TC*ARCL 

CQPF0N/C3/X0M(100I,Y0M( 100) ,S( 100 ) f SOL 1 100 ) » AE { 100 ) f TSE( 100) , 
1TAWL<100),TAWT(100),T8AR(100) tRW(lOO) tSWClCO) tSUTHL(lOO) , 
2RHSWU00) fRHSEt 100) ,HEADVil( 100 ) , HEADE ( 100 ) ,NUW( 100 ) t MUB AR ( 100 ) , 

3AA (ICO) ,68(100) tFF( 100) , DUDS ( 100 ), OMDS ( 100 ) ,OMDL ( IOC ) 

REAL NUSZ,NUSZ,MUTZ,NUTZ,ME,NUW,MUBAR 
INF =« 0 

81* l.+.E*(GAM-l.)*ME(l)**2 

82* l.+(2.79-1.78*PR**.5)»{ ( l.+SW( 1 ) )*81-1 . ) 

S3* -NUTZ*SUTHL( 1 )»ARCL/ATZ/DMOL ( I )*B1** ( ( 3.-GAM ) / ( 2. *GAM-2. ) ) 

BA* -1,1138*82 

85* 2.38All*82+( A.65*PR**( 1. /3. )-3.65*PR**.5)*PR**.5*(Bl-l, ) 

FX* (B3*XX)**,5*(84*XX+B5) 

IF (XX.EC.O.) GO TO 10 

DFX* ,5*(B3*XX)**(-.5)*B3*(BA*XX+B5)+BA*(B3*XX) **,5 
RETURN 
!1C INF * 1 

DFX * I.EIO 

RETURN 

END 
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REAL FUNCTION INTUXX) 

CONHCN/Cl/GAMjR, PTZ t TTZ , UPMACH, NST ,NVP ,NTURB , KP VM,KEM,KSMTH , 
ikspln»kle,katch,cthet,dlam,tlam,dturb,tturb,kpre,kgrao,ksoe,klam, 
2KMAINfKPPQF,X(100) ,Y( 100 ) ,PRES ( 100 J ,UE ( 100 ) , ME ( 100 ) ,POPTZ ( 100 » , 
3VCVCR(10C } ,TWALJ 100) 

C0PKCN/C3/X0M(100) ,Y0Ml 100 ) , S( 100) » SOU 100 ) f AE ( 100 ) t TSE ( 100 ) . 
ITAWLl 100),TAWT(100),TBAR( 100),RW(100) »SW(100) ,SUTHL(100) , 
2RHSW(100),RHSE( 100 ) i HEADWl 100 1 ,HEADE (100 ) ,NUW( 100) .MUBARdOO), 
3AA( ICO) ,EBUOO), FF( 100) , DUDS (100)* DMOSi 100 ) lOMDL ( 100 ) 
C0(^MGN/C5/SHAPL< 100 ) , SHAPK ( ICO ) * B,NS 
REAL ME,NUW,MUBAR,INT1 
CALL LGRNGE<SOL,ME,NST,XX,ANS) 

INT1= ANS**(B-1. )/( { l.+.5*(GAM-l.)*ANS**2)** 

1((3.»GAM-1, )/ (2.*GAM-2. ) ) ) 

RETURN 

END 


$iSFTC IMG2 DECK 

REAL FUNCTION INT2(XX) 

CCPf'CN/Cl/GAM,R, PTZ, TTZ, UPMACH, NST, N VP ,NTURB,KPVM,KeM,KSMTH, 
1KSPLN,KLE,KATCH,CTHET,0LAM,TLAM,DTUR8,TTURB,KPRE,KGRAD,KSDE,KLAM, 
2KMAIN,KPROF,X(100),Y(100 ),PRES( 100) ,UE(100),ME( 100) ,P0PTZ(100) , 
3V0VCRU0C ) ,TWAL( 100) 

C0MMCN/C3/X0M( 100) , YOM( 100) ,S( 100) , SOL ( 100 ) , AE ( 1 00 ) ,TSE(100) , 
1TAWL( 100),TAWT( 100) ,TBAR( 100),RW( 100) ,SW(100) ,SUTHL(100) , 

2RHSVV ( 100) ,RHSE( 100),HEA0W( 100 ) , HEAOE { 100 ) , NUW ( 100 ) ,MUBAR( 100) , 
3AA(1C0) ,681 100), FF( 100) , DUDS ( 100), DMDS( 100), DMDL( 100) 
CGNMCN/C5/SHAPL( 100 ), SHAPK ( ICO ), B , NS 
REAL ME,NUW,MUBAR, INT2 
IF (NS.LT.A) GO TO 10 
CALL LGRNGE(SOL, SHAPK,NS,XX, INT2) 

RETURN 

10 DO 20 J=2,NS 

IF (SCL( J ) .LT.XX ) GO TO 20 

INT2=? SHAPK(J-l)+( SHAPK! J)-SHAPK(J-l) )*(XX-S0L(J-1) )/( SOLI J)-SOL(J 
1 - 1 )) 

RETURN 
i2C CONTINUE 
RETURN 
END 
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SlfffTO SPLIN DECK 


SUERCUTINE SPLIN E(X,Y,-N,DYDX,02YDX2) 

SPLINE FITS A SPLINE CURVE TO X AND Y 

AND CALCULATES FIRST AND SECOND DERIVATIVES AT THE SPLINE POINTS 
END POINT SECOND DERIVATIVES EQUAL THOSE AT ADJACENT POINTS 

01 PENSION X(N),Y(N) ,DYDX(N)#D2YDX2IN) 

DIKENSION GIIOOI ,H(100) 

Gll)=? -1. 

M(l)» 0. 

Nl* N-1 

IR <M.LT,2) GO TO 20 
to 1C 1=2, N1 

AS {xm-x(i-in/6. 

(XtI+iy-X( I) )/6. 

C* 2.*( A+B)-A*G< I-n 

Q= (Yn+i)-Y( I) )/(x( i+i)-x( I) )-(Ym-Y( i-i)}/(x( n-xn-in 
6-m= B/C 

1C HI I )= ( C-A*H( I-l ) )/C 
2C D2YDX2(N»= H { N1 ) / ( 1. +G( N1 ) ) 

00 3C 1=2, N 
Ks N41-I 

30 D2YCX2(K)= H( K )-GI K ) *D2YDX2{ K+l ) 

DVDXm= (Xm-X(2) )/6.*(2.*D2YDX2(l)+D2YDX2(2) ) + (YI2)-Y(l))/(XI2) 

i-xnn 

DO 4C 1=2, N 

AO DYCXm= tx( n-Xn-l))/6.*(2.*D2YDX2n )+D2YDX2(I-l) ) + IYm-Y(I-l) ) 
i/(X( n-x< i-m 
RETURN 
ENC 
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$IBFTC RCC 


DECK 


SUBROUTINE ROOT! A»^B, Y.FUNCT, TOLERY,X,DFX ) 

ROOT FINOS t ROOT FOR (FUNCT-YI IN THE INTERVAL {A,B) 

XI* A 
X2* B 

GALL RUNCTIX1»FX1,DFX,INF» 

1C OQ 30 1*1,20 
X* (X1+X21/2. 

CALL FUNCT(X,FX,DFX, INF) 

IF ( (FXI-Y)*(FX-Y1.GT,0. ) GO TO 20 
X2= X 
GO TC 30 
20 XI* X 
fXl* FX 
30 CONTINUE 

IF (ABS(Y-FX) .LT.TOLERY) RETURN 
MRITE(6, ICQO) A,B,Y 
STCP 

lOOC FORMAT! ////4X,49HR00T HAS FAILED TO CONVERGE IN THE GIVEN INTERVAL 
1/4X,3HA *,G14.6, 10X,3HB =, G14.6 , lOX, 3HY *,G14.6) 

END 
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$IBFTG LGRNG DECK 


SUBRCUTINE LGRMGEPXiY^N, ARG, ANS) 

LGRNGE PERFCRMS 4 POINT LAGRAN6IAN INTERPOLATION 
01 PENSION X(N),YIN) ,XX(4),YY(4) 



IF 

{ARG-Xt2)) 10,10,20 

ro 

f¥K 

* 1 


00 

TO 70 

20 

IF 

(ARG-XI.N-D) 40,40,30 

30 

PM 

N-3 


GO 

TO 70 

AC 

N1 

= N-1 


DO 

60 1*2, N1 


IF 

(ARG-Xlin 50,50,60 

50 

m 

* 1-2 


GO 

TO 70 

^0 

CONTINUE 

TO 

DO 

80 1=1, A 


PPP 

' * MM+I-1 


XXII) = XIPNM) 

80 

YYII) = YIPMM) 


Cl => l(ARG-XXr21 )*IARG-XX(3))*<AR6-XX{4)))/ 
l!XX(U-XXl2n/(XX( 1)-XX(3M/(XX( 1»-XX(4) ) 

C2 = ((ARG-XXmH‘(ARG-XX(3n*(ARG-XX(4)))/ 
1CXX(2)-XX( 1) )/(XX<2)-XX(3n/(XX(2I-XX(4)) 

C3 » K ARG-XXlin*(ARG-XX(2»»*(ARG-XX(4) n/ 
1IXX(3I-XX( n)/(XX(3)-XXI2) )/(XX(3)-XX(4) ) 

C4 * U ARG-XXn) l♦(AR6-XX<2n♦(ARG-XX^3) ))/ 
ltXX(4l-XXn»)/(XX{4)-XX(2n/(XX(4)-XX(3) ) 
ANS a C1’*YY(1)+C2=»YY(2)+C3*YY(3)+C4*YY(4) 
RETURN 
END 
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$1BFTG CURVF DECK 


SUBRCUTINE CURVFKCOEF, ANS,X,Y,NX,NY) 

EVALUATE THE POLYNOMIAL FUNCTION, ANS*F(X,Y), USING COEFFICIENTS, COEF 

BIMENSION CQEF(20) 

NXl = NX41 
NYl * NY41 
ANS =f CCEF(l) 

If (X.EQ..O.AND. Y.EQ..O) RETURN 
IF <Y.EQ..C) GO TO 10 
IF <X.EC..O) GO TO 30 
GO TC 50 
10 DG 20 1=2, NXl 
20 ANS = ANS4C0EFI I )*X**n~l) 

RETURN 

30 DO 40 I=2,NY1 
K ~ (I-D4NX1+1 
40 ANS = ANS4C0EF(K l*Y**t I-l) 

RETURN 
50 ANS =* .0 

DO 60 1=1, NYl 
DC 60 J=1,NX1 
K = ( I-1)«NX1+J 

60 ANS =; ANS+C0EF(K)*Y**(I-1)*X**C J-1) 

RETURN 

END 


SIBFTC SIMP DECK 

FUNCTION SIMPS1CXI,X2,FUNC,KSI6) 

DIMENSION V(200) ,H ( 200) , A{ 200 ) , BI 200 ) ,C ( 200 ) , P ( 200 ) ,E ( 200 ) 

LOGICAL SPILL 

DOUBLE PRECISION ANS,Q 

DATA TWO, THREE, FOUR, THIRTY/2. 0,3. 0,4. 0,30.0/ 

DATA T,NMAX,NSIG/3.0E-5,200, 1/ 

C INITIALIZE FIRST ELEMENTS OF ARRAYS. 

V=X1 

H=(X2-V)/TW0 
A=^FUNC( V) 

B=FUNC( V4H) 

CfcFUNCt X2) 

P=H*(A+FCUR*B+C) 

E=P 

ANS=P 

N=1 

FRAC=T 

SPILL=5. FALSE. 

1C TEST=ABS(FRAC*ANS) 

K'=N 

DO 30 I=1,K 

C TEST MAGNITUDE OF 4TH ORDER ERROR IN THIS INTERVAL. 

IF (ABS(E( m.LE.TEST) GO TO 30 
IF (N.LT.NMAX) GO TO 20 
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C GC TO FINISH IF STORAGE IS FILLED UP. 

SPILL=;.TPUE. 

KSIC=KSIGHNSIG 
GC TO 40 

C 'SLBDIVIEE INTERVAL AGAIN TO REDUCE 4TH ORDER ERROR. 

2C N=N+1 

VlNHVl IHH( I) 

MTN) = H( I I/TWO 
A(N):*B( I ) 

B(M=»fiUNC(V{N»+H{N) ) 

C(M*C( I » 

P<N)»H(N)*IA(N)+FOUR*B{N J+CINH 
HlilMHIN) 

Btn*FUNciv(n+Hm) 

crn^A(N) 

Q*PC IJ 

ptiHHi n*iA( n+FouR*B< n+c( in 

Q=Pl mP(N)-Q 
ANS=ANS+C 
El I)-C 
ELN)=<J 
3C CONTINUE 

C TEST ALL INTERVALS AGAIN IF ANY WERE SUBDIVIDED THE LAST TIME, 
IF (N.GT.K) GO TO 10 
>40 QnC.C 

DO 5C 1=1, N 

50 Q^c+em 

C TIGHTEN ERROR LIMIT IF TOTAL ACCUMULATED ERROR TOO LARGE. 

IF (ABS(C/T).LE.A8S(ANSI. OR. SPILL! GO TO 60 

FRAC=FRAC/TWO 

GO TC 10 

C FINISH CALCULATION. 

'60 SIMPS1=( ANS + Q/THIRTY)/THREE 
RETURN 
ENC 
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SIB6FTC AFTf»IX DECK 


SUBROUTINE AFMIX I ALPHl, TEf SP t XMFSl ,NP J 
DIMENSION XXXIIOO) ,YYY( 100) i SSSI 1001 

CGMMCN/Cl/GAM,RtPTZ,TTZ,UPMACH»NST,NVP,NTURB,KPVM,KEM,KSMTHt 
1KSPLN,KLE,KATCH,CTHET,DLAM,TLAM,DTUR8,TTURB,KPRE,KGRAD,KSDE,KLAM, 
2KMAIN,KPR0F,X(100),Y( 100 ) ,PRES ( 100) ,UE ( 100 ) t ME! 100) tPOPTZIlOO) , 
3YGVCR(1CC) ,TWAL< 100) 

GOMMON/C3/XON(100),YOMI 100 ) t SI 100) t SOL 1 100) ,AE(100) ,TSE(100) , 
ITAWLIlOOIfTAWTI 100),TBAR( 100),RW(100),SW(1CO) fSUTHLIlOO) , 

2RHSWI 100) ,RHSE( 100) ,HEAOW( 100) tHEADEI 100) ,NUWI 100) ,MUBAR ( 100 ) , 
3AAI100) fBB( 100) f FF( 100), DUDS (100) »OMDS( 100) tOMOLI ICO) 
CCMMCN/CA/THETI 100), DEL SRI 100), DELTA 1 100), FORM (100) , 

IFORMItlOC ) ,FORMTRI 100),RTH( 100),RTHII100),CF( 100), 

ITAUWI 100 ) ,NUSSI 100) ,DTDY( 100 ) ,HTRANI 100 ) ,CRN ( 100 ) 

ECU I VALENCE I X, XXO{ ) , ( Y, YYY ) , I S, SSS ) , INST,N ) 

REAL ME 

KGCE = 0 

0BLS * DELSRIN) 

DELP a DELSRINP) 

THETS = THETIN) 

THETP = IHETINP) 

WRITEt6,5) 

5 FORMAT! in, 20X,22HAFTERMIXING PROPERTIES) 

MRITEI6,n 

1 FORMAT! 11-0, 10X,40HN0ZZLE WITH NO BOUNDARY LAYER CORRECTION) 

2 VVCR1=^SQRTI I (GAM-H.0)/2.0*XMFS1**2)/I1.0*(6AM-I.0)/2,0*XMFS1**2)) 
XX*SP*COSI ALPMl) 

OELSRTa (DELS*DELP)/XX 
THETA*^ (THETS*THETP)/XX 
DTE a TE/XX 

AFSl a (6AM-1.0) /IGAM+1.0)*VVCRl**2 
A = l.O-CELSRT-OTE-THETA 
Al= l.O-CELSRT-OTE 
IF (A .LE. 0.0) GO TO 16 

C=iI(1.0-AFSl)*(GAM + 1.0)/(2.0*GAM) + ICOS(ALPHl) )**2* 

1A*VVCR1**2 )/ICOS(ALPHl)*Al+VVCRl) 

D = VVCR1#SINIALPH1)*A/A1 

VXVCR2 = GAM*C/ (GAM*1.0)-SQRTI (GAM*C/(GAM+1.0) )**2 
l-l.C+:iGAM-1.0)/l GAM + 1.0)*D**2) 

@0 TC 4 

3 KGCEa2 
«RITEt6,10) 

10 FORMAT! 1FO,10X,19HSUPERSONIC SOLUTION) 

VXVCR2 * GAM*C/(GAM+1,0)+SQRTI IGAM*C/(GAM+1.0) )**2 
l-1.0+6GAM-1.0)/( GAM+1,0)*D**2) 

4 BENCH = (GAM-1. 0)/(GAM+1.0)*ID**2+VXVCR2**2) 

IF (CENCF.GE. 1.0) RETURN 

0ENSR2= 1 1.0- (GAM-1.0 )/(GAM+1.0)*ID**2+VXVCR2**2)) 
1*4(1.0/(GAM-1.0) ) 

BENSRl a 1.0/((1.0+(GAM-1.0)/2.0*XMFSl**2)**(1.0/( 

I6AM-1.0) )) 

PR2 a DENSR2**GAM 

PT2PT0=(CENSR1*VVCR1*C0S(ALPH1)*A1)/(0ENSR2*VXVCR2) 

PTCP2 a 1 .0/(PT2PT0*PR2) 

EBAR2 a ((PT2PTO)**((1.0-GAM)/GAM)-1.0)/IPTOP2** 
l«tGAM-l.C)/GAM)-1.0) 

ETAN a 1.0-EBAR2 

VVCR2 a SQRT(D**2+VXVCR2**2) 

XM2 a SQPTI I ( 2.0/1 GAM+1.0) )*VVCR2**2)/( 1.0- ( (GAM-1. 0)/(GAM+1.0) ) 
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I+VVCR2**2)) 

T2TTC = 1.0 - <GAW-1.1/(GAM+1.)*VVCR2**2 
ALPH2 = ATAN<0/VXVCR21 
XRXl =» XPFS1*C0S(ALPH1) 

XPX2 XP2*C0S< ALPH2) 

ALPH2 * ALPH2*57.2958 
ALPHl = ALPH1*57.2958 

WRITE (6,6) XMFSl ,SP,TE,XM2,VVCR1,XMX1,XMX2 
5 FGRPATt U-0,7HXMFS1 =,F6.4,2X,9HSPACING *,F7.5,2X,4HTE *,F7.5,2X, 

15HXK2 =,F6.4,2X,8HV/VCR1 =,F6.3,2X,6HXMXl = ,F6. 3 ,2X ,6HXMX2 =,F6.3) 
WRITE(6,7)ALPH1, ALPH2,PT2PT0,PT0P2,T2TT0,VVCR2,EBAR2,ETAN 
7 FGRPATC 1F0,6HALPHI*,F7.3,2X,6HALPH2=,F7.3,2X,8HPT2/PT0=,F7.4,2X, 

17HPT0/P2»»F9.3,2X,7HT2/TT0=,F7.4,2X,7HV/VCR2*,F6.3,2X,6HEBAR2*, 
2F7.5,2X,6HETA-N*,F6.4) 

IF (KCDE .EG. 2 ) RETURN 

IF (KCDE. EC. 1 ) WRITE(6,15) DELSRI N+1 ) ,DELSR ( N + 2 ) , THEKN+l ) , 
ITHET(N*2I 

15 FCRPATt 1F0,12HDELSR(N*1) F8.6, 2X , 12HOELSR (N+2 ) =,F8.6,2X, 

111HTHET(A+1) «,F8.6,2X, llHTHET(N+2) *,F8.6) 

IF (KQDE .EG. 1 .AND. XMXl .LT. 1.0000 ) RETURN 
ALPH = AtPHl*. 017453 
IF( KCDE. EC. 1 ) GO TO 3 

16 SP= 2.0*(VYY(NP)+DELP)/COSCALPHl) 

CISPLACEMENT THICKNESS HAS BEEN ADDED TO THE NOZZLE EXIT 

INCREASING THE STRAIGHT SECTION 
Z » 2.0*(YYY(NP)+DELP)*TAN(ALPH1) 

CELZ Z-(XXX(N)-XXX(NP) ) 

XXX(N + 1) =» CELZ + XXX(N) 

¥YY(N+1) * YYY(N) 

FURTHER INCREASING OF THE STRAIGHT SECTION DUE TO INCREASING 
BOUNDARY LAYER 

NOZZLE ANCLE IS CONSTANT 

FA » ATAN((DELSR(N)-DELSR(N-2))/(SSS(N)-SSS(N-2) )) 

OELSRIN+l) * DELSR(N)+DELZ*TAN(FA) 

FAPT = ATAN((THET(N)-THET(N-2) )/(SSS(N)-SSS(N-2) ) ) 

THET(N+1) * THET (N)+DELZ*TAN(FAMT) 

ALPVI s 3.14159 - ALPHl - (3.14159/2.0 +FA) 

DS * SIN(3. 14159/2. 0+FA)*(DELSR(N+l)-DELSR(NP) ) /SIN(ALPW) 

BF = SIN(ALPHl)*(OELSR(N*l)-DELSR(NPn /SIN(ALPW) 

0ELZF = CF*COS(FA) 

OELSRlN+2) = DELSR(N+1) + DELZF*TAN ( FA ) 

THET(N+2) * THET(N+1) + DELZF*TAN( F AMT ) 

XXX(N + 2) XXX(N + 1) + DELZF 

YYY(N*2) YYY(N) 

SP ^ SP 4DS 
DELS - DELSR(N+2) 

THETS = THET(N+2) 

WRITE(6,20) 

20 F'CRMAT( 1H0,10X,37HN0ZZLE WITH BOUNDARY LAYER CORRECTION) 

KGCE =! 1 
GO TC 2 
END 
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SIBFTC HCI2C DECK 


SUBROUTINE NOZZLCINSPtALPHl.NCP ) 

DIPENSION YYCU00),YYCL(100) ,XXT( 100 »,YYT( 1001 tXXTU 100) ,YYTH 100) 
DIMENSION XXX(100),YYY( 100) 

eOPPCN/Cl/GAM,R,PTZtTTZ,UPMACM,NST,NVP,NTURB,KPVM,KEM,KSMTH, 
lKSPLN»KLE,KATCHtCTHET,DLAM,TLAH,DTUR8,TTURB,KPRE,KGRAD,KSDE,KLAM, 
2KPAIN»KPROF,XUOOJ,Y( 100),PRES( 100 ) ,UE ( 100 ) ,ME 1 100 ) ,POPTZ 1 100) , 
3V0VCRC10C ) ,TWAU 100) 

CQPPCN/CA/THET(100),DELSR( 100) t DELTA! 100) t FORM! 100 ) t 
liORPH10C),FORMTR(100),RTH(100) »RTH1(100),CF!100), 

ITAUWI 100),NUSS(100),DTDY(100),HTRAN!100),CRNI100) 
soul VALENCE ( X, XXX ) , ( Y, YYY ) ♦ (NST ,N ) 

Ni ^ N+2 
X3.UB ^ XXXINSP) 

BO 1C I » NSPrNl 
XXX! I ) =XXXn ) -XSUB 
¥YC( I) « YYY! n + DELSR! I ) 

IF ! I.GT.NCP) GO TO 10 
YYCL(I) * -YYC! I ) 

10 CONTINUE 

ALPHl * ALPH1+. 017453 
DO 11 I=NSP,N1 

XXT! II * XXX! I)*COS( ALPHD- YYC ! I ) *SIN ! ALPHl ) 

YYT! I ) ■ YYC! I)*COS! ALPHl) +XXX! I )*SIN! ALPHl ) 

IF ! I.GT.NCP) GO TO 11 

XXTL! I)= XXX! I)*C0S1 ALPHD-YYCL! I )*SIN!ALPH1) 

YYTLII)»YYCL( I) *COS( ALPHl) 4XXX! I )*SIN! ALPHl ) 

11 CONTINUE 
WRITE!6,1) 

1 FGRPAT! IH,20X,18HN0ZZLE COORDINATES, IX, 18HSUPERS0NIC SECTION) 

NRITE(6,2) 

2 FORMAT! IFO , lOX, 37HN0ZZLE WITH BOUNDARY LAYER CORRECTION) 

WRITE!6,3) 

3 FORMAT! IFO , 9X , IHX , 1 IX , IHY, lOX , 3HXTU, 9X , 3HYTU ,8X , 3HXTL , lOX ,3HYTL ) 

DO 12 I=NSP,N1 

12 MR ITS C6, 4 ) XXX! I ) ,YYC1 I), XXT! I), YYT! I ) , XXTL ! I ) , YYTL 1 1 ) 

4 FORMAT! IF ,6F12.6) 

RETURN 

ENC 


Lewis Research Center, 

National Aeronautics and Space Administration, 
Cleveland, Ohio, May 19, 1971, 

120-34. 
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APPENDIX A 


PROGRAM CHANGES FOR A GAS OTHER THAN AIR 

The program gas properties are set up for air. This program can be easily 
changed so that it applies to gases other than air . The changes are all made in sub- 
routine PRECAL, 

The coefficients read in by the DATA statements for CMU, CPR, and CTC arrays 
must be changed. The equations for the curve fits have the following form: 



If the munber of coefficients changes from five in any case, this must be reflected both 
in the DATA statements and later in the calls on CURVFT where CMU, CPR and CTC are 
used. If the properties are put in a different form, these cards must be removed. 

The sea- level reference values in U.S. Customary and SI Units of temperature 
(TSLE, TSLM), viscosity (MUSLE, MUSLM), and thermal conductivity (TCSLE, 

TCSLiM) must be changed. 

The value of Sutherland's constant (TCON) and the computation of (SUTHL(I)) 
will have to be changed. A temperature-viscosity law of the following form was used; 



where for air 


ksu=l 


^rjM/2 /t^ + TCON^ 


T + TCON 


with TCON being Sutherland's constant for air. 
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APPENDIX B 


ADDITIONAL OUTPUT 

This output is obtained when KPRE, KGRAD, KSDE, KLAM, and KPROF are set 
equal to 1. 

The output corresponding to KPRE are the geometric variables 
X X” coordinate 5 m (ft) 

Y Y- coordinate, m (ft) 

S surface length, x, m (ft) 

XOM ratio of X to maximum X-coordinate 
YOM ratio of Y to maximum X-coordinate 
SOL ratio of surface length to total arc length 

The next part of this output gives the local speed of sound and the several temperatures 
AE local free -stream speed of sound, m/sec (ft/sec) 

TSE static temperature, K (°R) 

TWAL wall temperature, K (°R) 

TAWL laminar recovery temperature, K (°R) 

TAWT txirbulent recovery temperature, K (*^R) 

TEAR reference temperature, K (°R) 

The final part of this output gives 

RW Reynolds number at wall, RW = (UE)(S)/NUW 

SW temperature function at wall 

SUTHL value of coefficient in Sutherland's viscosity temperature formula 

O O 

RHSW static density based on wall temperature, kg/m (slug/fr) 

O O 

RUSE static density based on free-stream temperature, kg/m (slug/ft ) 

HEADW velocity head based on density at wall, N/m^ (Ibf/ft^) 

HEADE velocity head based on free-stream density, N/m^ (Ibf/ft^) 

p p 

NUW kinematic viscosity at wall, m /sec (ft /sec) 

p p 

MUBAR dynamic viscosity based on reference temperature, (N)(sec)/m ((lbf)(sec)/ft ) 
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Output corresponding to KGRAD contains the three gradients of velocity and Mach 
number along the surface computed by the spline curve-fit method: 

DUDS dUE/dx, sec"^ 

DMDS dME/dx, m"^ (ft“^) 

DMDL dME/dSOL 


Output corresponding to KSDE contains the niunerical solution of the laminar and 
turbulent differential boundary- layer equations. In the laminar case, the solution is the 
correlation munber (CORLN). In the turbulent case, the solution is the incompressible 
form factor (FORMI), and a function (F) of the momentum thickness. These solutions 
are printed with respect to the surface length (S). 

Output corresponding to KLAM contains the variables used in the laminar subrou- 
tine to check for the position of instability and transition. RTHI, increasing from sta- 
tion to station, and RCRIT and RTRAN, decreasing from station to station, are used in 
this analysis. When RTHI grows larger than RCRIT, instability has occurred. When 
RTHI bypasses RTRAN, transition is assvuned to occur. For most cases of supersonic 
nozzle flow, the flow is laminar. The variables listed are 


CORLN correlation number 

SHEAR shear parameter 

DTH ratio of transformed displacement thickness to transformed momentum 

thickness 

FORMTR transformed form factor 

SHAPE Pohlhausen shape factor based on boundary- layer thickness 
RTHI incompressible momentum -thickness Reynolds number 

SHAPK dimensionless shape factor based on momentum thickness 
RCRIT . critical incompressible momentum -thickness Reynolds number 
KBAR mean shape factor based on momentum thickness 

DIFF difference between transition and instability momentum-thickness Reynolds 

numbers 

RTRAN incompressible momentum-thickness Reynolds number used in checking for 
transition point 

Output corresponding to KPROF contains the velocity profiles at each station along 
the surface. The output listed is 

Y/DELTA ratio of distance normal to sirrface in y-direction in boundary- layer profile 
to boundary-layer thickness 
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Y 

y/XMAX 

u 

U/UE 


distance normal to surface in y-direction in boundary- layer profile, m (ft) 

ratio of Y to maximum X- coordinate 

velocity within boundary layer, m/sec (ft/sec) 

ratio of U to free-stream velocity 
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